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1.  INTRODUCTION 


Regenerative  liquid  propellant  gun  (RLPG)  technology  is  being  considered  in  the 
development  of  the  next  U.S.  Army  artillery  weapon.  Test  fixtures  in  30-mm,  105-mm,  and 
155*mm  sizes  have  been  built  and  fired.  The  data  from  ail  these  fixtures  have  been 
extensively  analyzed  to  better  understand  the  RLPG  process  (Coffee,  Wren,  and  Morrison 
1989, 1990;  Wren,  Coffee,  and  Morrison  1990;  Coffee  et  al.  1991).  Previous  reports  describe 
the  design  of  the  Concept  VIC  gun  fixtures  and  the  modeling  modifications  pertinent  to  this 
design.  Agreement  between  the  interior  ballistics  simulations  and  the  experimental  data  is 
quite  good.  Details  of  the  lumped  parameter  gun  code  used  for  these  past  simulations  have 
been  described  in  previous  publications  (Coffee  1985,  1988;  Morrison  and  Coffee  1990). 

A  Concept  VIC  RLPG  is  shown  in  Figure  1.  An  external  solid  or  liquid  propellant  igniter 
venting  into  the  combustion  chamber  initiates  the  ballistic  cycle.  The  chamber  pressure  forces 
both  the  control  and  injection  pistons  rearward.  Liquid  propellant  is  then  injected  from  the 
liquid  reservoir  through  the  annulus  between  the  pistons  into  the  combustion  chamber  where  it 
bums.  The  motion  of  the  control  piston  Is  modulated  by  the  damper  assembly. 

Almost  all  firings  of  liquid  propellant  guns  have  shown  high  frequency  oscillations.  As  an 
example.  Figure  2  shows  pressures  from  shot  65  of  a  first  generation  155-mm  Concept  VIC 
regenerative  liquid  propellant  gun.  The  top  figure  shows  the  pressure  measured  at  the 
combustion  chamber  wall.  The  bottom  figure  (moved  down  100  MPa)  shows  the  pressure  at 
a  gauge  1 .7  m  downbore.  Note  that  large  pressure  oscillations  persist  for  a  long  distance 
down  the  gun  tube. 

Figure  3  shows  the  Fourier  transforms  of  the  two  pressures  over  a  time  window  from  14  to 
15  ms.  One  would  normally  expect  the  frequencies  to  be  near  the  natural  acoustic 
frequencies  of  the  chamber.  However,  for  such  a  large  chamber,  the  acoustic  frequencies  are 
small  (first  radial  ■  4.0  kHz).  It  is  not  dear  what  causes  the  high  frequency  oscillations. 

The  barrel  gauge  shows  frequendes  similar  to  the  chamber  gauge,  although  the  higher 
frequendes  are  missing.  Pressure  oscillations  are  driven  by  the  energy  released  by  the 
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Figure  3.  Fourier  Transform  of  the  Experimental  Data  From  14  to  15  ms.  E222  (line). 

B3  (dot). 

combustion  process.  Since  most  of  the  combustion  is  expected  to  occur  in  the  chamber 
(based  on  analysis  of  experimental  data),  it  is  also  unclear  why  there  are  still  large  magnitude 
pressure  oscillations  so  far  downstream. 

2.  MODEL  ASSUMPTIONS 

To  study  these  problems,  a  two-dimensional  computational  fluid  dynamic  model  of  the 
combustion  chamber/gun  tube  has  been  created.  At  the  present,  no  attempt  has  been  made 
to  model  the  liquid  reservoir  or  the  damper. 

In  the  Concept  VIC  regenerative  liquid,  the  propellant  is  initially  behind  two  pistons.  The 
inner  or  control  piston  moves  first,  opening  up  an  annular  vent  The  outer  or  injection  piston 
then  slightly  trails  the  inner  piston.  Propellant  is  injected  into  the  combustion  chamber.  Most 
of  the  combustion  takes  place  in  the  combustion  chamber.  The  gas  then  flows  into  the  gun 
tube.  There  is  a  large  area  change  from  the  chamber  to  the  tube.  In  the  30-mm  and  105-mm 
fixtures,  there  is  a  sharp  comer.  In  the  155-mm  fixtures,  there  is  a  taper  from  the  chamber  to 
the  gun  tube. 
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In  the  two-dimensional  model,  the  combustion  chamber  is  idealized  as  a  cylinder.  The  left 
wall  is  allowed  to  move  in  an  arbitrary  fashion,  roughly  mimicking  the  motion  of  the  pistons 
that  inject  liquid  into  the  chamber.  The  grid  is  attached  to  the  left  wall  and  stretches  as  the 
wall  moves.  The  individual  grid  points  are  moved  so  as  to  keep  the  grid  uniform  in  the  axial 
direction.  There  is  a  circular  hole  in  the  center  of  the  right-hand  wall  corresponding  to  the 
entrance  to  the  gun  tube. 

The  injection  process  is  simplified.  Liquid  droplets  of  a  specified  size  are  injected  at  a 
specified  rate  through  an  annulus  in  the  left  wall.  The  size  of  the  annulus  may  vary  with  time. 
Infinite  drag  is  assumed,  so  the  gas  and  liquid  at  a  given  point  have  the  same  velocity,  and 
only  one  set  of  momentum  equations  is  required.  The  droplets  then  burn  and  release  hot  gas, 
according  to  a  pressure  dependent  burning  rate. 

The  chamber  may  be  run  in  isolation,  with  exit  conditions  specified  at  the  gun  tube 
entrance.  However,  a  gun  tube  model  is  available.  The  gun  tube  is  also  considered  to  be  a 
cylinder,  with  the  projectile  at  the  right  end.  The  standard  equations  of  motion  for  the 
projectile  are  implemented.  As  the  projectile  moves,  the  grid  stretches.  Since  the  gun  tube 
length  increases  dramatically  during  the  course  of  the  firing  cycle,  additional  grid  points  are 
also  added  as  the  projectile  moves.  Values  on  the  new  grid  are  found  by  interpolation. 

One  option  allows  the  use  of  the  output  from  the  lumped  parameter  gun  code  (Coffee 
1985,  1988)  as  an  input  into  the  two-dimensional  code.  That  is,  the  velocity  of  the  left  wall 
(piston  velocity),  the  injector  area,  the  propellant  injection  rate,  and  the  injected  droplet  size 
are  ail  read  in  as  a  function  of  time  from  the  lumped  parameter  code  output.  This  allows  the 
two-dimensional  code  to  mimic  a  gun  firing  cycle.  In  addition,  a  primer  model  is  required. 

This  causes  hot  gas  to  be  injected  through  an  annulus  in  the  right-hand  wail  of  the  chamber  at 
a  specified  rate  at  the  start  of  the  integration. 

3.  NUMERICAL  PROCEDURE 

An  orthogonal  grid  is  set  up  in  the  combustion  chamber.  This  divides  the  chamber  into 
annular  control  volumes.  The  scalar  quantities  (pressure,  temperature,  and  density)  are 
assumed  to  be  uniform  within  a  control  volume.  The  axial  and  radial  velocities  are  defined  on 
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the  vertices  of  the  control  volumes.  The  position  of  the  vertices  may  be  arbitrarily  specified  as 
a  function  of  time.  This  type  of  mesh  is  called  an  arbitrary  Lagrangian-Eulerian  (ALE)  mesh 
(Hirt,  Amsden,  and  Cook  1974;  Cloutman  et  al.  1982;  Amsden  et  al.  1985)  and  is  useful  for 
representing  moving  boundaries.  In  this  case,  the  vertices  only  move  axially  with  the  left  wall, 
and  the  procedure  is  somewhat  simplified  from  a  general  ALE  formulation.  The  governing 
equations  are  cast  in  integral  form.  The  procedure  is  inherently  conservative.  That  is,  if 
mass,  momentum,  or  energy  leaves  one  control  volume  through  a  boundary,  it  automatically 
enters  the  neighboring  control  volume. 

An  explicit  numerical  integration  scheme  is  used.  Implicit  methods  improve  the  stability  of 
the  integration  and  allow  larger  time  steps  to  be  used.  However,  the  purpose  here  is  to 
actually  track  pressure  waves.  This  implies  that  very  small  time  steps  must  be  taken  to 
accurately  resolve  the  motion  of  the  pressure  waves.  Therefore,  the  additional  overhead  of  an 
implicit  method  is  not  useful. 

Turbulent  transport  processes  are  described  by  a  k-e  turbulence  model  (Jones  and 
Launder  1972;  Jones  and  Whitelaw  1982).  This  introduces  two  new  partial  differential 
equations  for  the  turbulent  kinetic  energy  and  the  dissipation  rate.  The  turbulent  viscosity  is 
computed  from  these  quantities.  There  are  much  simpler  turbulence  models.  However,  the 
k-e  model  is  the  simplest  model  that  allows  for  transient  effects.  That  is,  the  turbulent 
viscosity  can  start  out  negligibly  small  and  grow  during  the  firing  cycle.  More  complicated 
Reynold's  stress  models  are  judged  not  to  be  worth  the  additional  effect.  All  of  the  turbulence 
models  involve  empirical  parameters  that  must  be  set  by  comparison  with  experiment.  Since 
these  parameters  have  been  chosen  based  on  incompressible  flow  experiments  near 
atmospheric  pressure,  extrapolating  to  gun  conditions  is  only  an  approximation  in  any  case. 

Physically,  the  gas  velocity  at  the  wall  should  exactly  match  the  wall  velocity  (no  slip 
condition).  The  flow  very  near  to  the  wall  is  laminar.  Further  away  from  the  wall,  the  flow  will 
become  turbulent.  But  there  will  be  a  boundary  layer  where  the  flow  is  primarily  parallel  to  the 
wall.  Particularly  for  the  large  pressure  gradients  in  a  gun  simulation,  boundary  layers  at  the 
walls  will  be  very  thin.  It  is  impractical  to  use  enough  grid  points  to  actually  resolve  the 
boundary  layers.  Boundary  layer  theory  (law  of  the  wall)  is  used  to  set  the  boundary 
conditions  (Bradshaw  1978).  That  is,  the  expected  velocity  at  the  edge  of  the  boundary  layer 
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is  used  as  the  wall  velocity  in  the  code.  Much  of  the  turbulence  is  generated  near  the  wall, 
and  this  can  also  be  computed  using  boundary  layer  theory. 

3.1  Governing  Equations  -  Laminar  Flow.  The  governing  equations  are  given  in  Cartesian 
coordinates.  While  cylindrical  coordinates  might  seem  more  appropriate  for  the  problems  of 
interest,  it  is  more  convenient  to  consider  the  geometry  when  the  equations  are  recast  into 
finite  volume  form.  The  continuity  equation  for  two  dimensions  is  (Bird,  Stewart,  and  Lightfoot 
1960) 


3p/3f  +  3(p  vx)ldx  +  3(p  v^/Sy  -  0  .  (1) 

Now  consider  some  volume  in  three-dimensional  space.  By  the  divergence  theorem,  the 
above  equation  can  be  written  as 

d/dtfvpdV  -  -fvp(pvx)/dx  +  d(pvy)/dy]dV 
-  +  p  vyj)ndS. 


The  first  two  integrals  are  volume  integrals.  The  last  integral  is  a  surface  integral.  The 
vectors  /  and  /  are  the  unit  vectors  in  the  x  and  y  direction,  and  n  is  the  unit  outward  normal 
vector  to  the  surface.  So  any  change  in  the  mass  in  the  control  volume  is  due  to  mass  flux 
across  a  boundary. 

The  energy  equation  can  be  written  as  (Bird,  Stewart,  and  Lightfoot  1960) 

p  cy3773f  ■  -p  c,v,3T/3x  -  pcYvy  dT/dy 
-  T(3p/3T)p(3yr/3x  +  3v/3y). 


The  viscous  heating  term,  which  is  almost  always  small,  is  eliminated.  The  thermal 
conductivity  term  is  also  not  included.  For  the  problems  of  interest,  the  temperature  is  fairly 
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constant,  and  the  thermal  conductivity  term  is  much  smaller  than  the  convection  terms.  For 
the  Noble-Abel  equation,  the  derivative  of  p  with  respect  to  7  with  the  density  constant  is  just 
pressure  over  temperature.  Making  this  substitution  and  dividing  by  the  specific  heat, 

pdT/dt  »  -p  vgdTldx  -  p  vydT/dy 
-  p/Cv(dVK/dx  +  dvy/d y). 


Using  the  continuity  equation,  this  can  be  written  as 

3(p7)/a/«  — 3( p  v„  T)/dx  -  d(p  vyT)/dy 

-p!  cv(dvx/dx  +  dvy/d  y).  (2) 


The  divergence  theorem  is  applied.  The  pressure  is  assumed  to  be  constant  in  the  control 
volume  and  can  be  taken  outside  the  integral.  Then 


9/dtfjp  T)  d V  -  - /s(p  77  .  p  v, 77 )  ndS 
-  (p/c,)  J"s(v,/  ♦  vtJ)ndS. 


The  basic  variable  is  the  product  of  the  mass  and  the  temperature.  This  can,  as  before,  be 
changed  by  convection  (first  right-hand  term).  Hot  or  cold  gas  can  enter  through  a  boundary. 
The  temperature  can  also  be  changed  through  work  (second  right-hand  term).  Note  that  in 
the  last  term  the  pressure  is  not  evaluated  at  the  boundary  but  is  the  pressure  in  the  control 
volume. 

The  momentum  equation  in  the  axial  direction  is  (Bird,  Stewart,  and  Lightfoot  1 960) 

p  dvx/dt  *  -pv^dvjdx  -  pVyBvx/dy  -  dp/dx 
-  dxjdx  -  dxy,/dy. 
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Using  the  continuity  equation,  this  can  be  rewritten  as 


3(pv-,)/3f  «  -3(pv*)/3x  -  3(pvjrv,r)/3y  -  3p/3x 
-  3t„/3 x  -  3t/jr/3y. 

For  a  Newtonian  fluid,  the  shear  stresses  can  be  written  as 

x„  »  —  4/3  jx 3 1^/3 +  2J3\iBvyldy 
xyy  •  -VZiidvJdy  +  2/3p3v,/3x 
x,,  -  V  -  -n(3  v,/3y  +  dvjdx) , 

so 

3(P  v'Jr)/3^  ■  -  9(PKra)/9>f  -  3<p v,  v^J/ay  -  3p/3x 
+  3/3x(4/3p.3v,/3x  -  TJZyidvfly) 

+  3/3y(p.3v’jr/3y  +  pSv^x). 


Applying  the  divergence  theorem, 

3/3fJ^(p  v,)dV  -  -  Js(pv,y,/  +  p  v,^/  +  pl)ndS 
♦  JJ(4/3n3v,/3x  -  2/3n3v/3y)/ 
+  (p.3v/3y  +  p.3vy/3x)/]ndS. 
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Similarly,  the  momentum  equation  in  the  radial  direction  is 


p9  Vy/dt  »  -  P  V^Vy/dX  -  QVydVy/dy  “  ^  p!  8 
-  dxgy/dx  -  dXyyldy. 


Using  the  continuity  equation  and  replacing  the  shear  stresses, 

8(p vy)/dt  -  -  d(pvxvr)/dx  -  8(pvyvy)/8y  -  dp/dy 
+  3/8x(p.8v,/8y  +  \Ldvy/dx) 

+  8/8y(4/3p.8vy/8y  -  ZJZpdvJbx) .  (4) 


Applying  the  divergence  theorem. 


8/8 1  jv{  p  vy )  dV  -  -  Js  ( p  vM  vy  /  ♦  p  vy  vyJ  +  pj )  n  dS 
*  +  \idvy/dx)i 

+  (4/3^3 vy/8/  -  2J3\idvg/dx)  /]  ndS. 

Values  for  the  viscosity  of  the  fluid  are  required.  The  liquid  is  assumed  to  have  a  negligible 
effect  on  viscosity,  and  only  the  gas  viscosity  is  computed.  The  most  common  propellants 
used  are  LGP1845  and  LGP1846  (Decker  et  al.  1987).  For  simplicity,  the  assumption  is 
made  that  LGP184S  goes  completely  to  the  major  species  C02,  N2,  and  H20  (no  minor 
species).  The  viscosity  of  these  species  as  a  function  of  temperature  can  be  found  using 
Lennard-Jones  (Reid  and  Sherwood  1 966)  (non-polar)  or  Stockmayer  (Monchick  and  Mason 
1981)  (polar)  parameters.  For  the  mixture  viscosity,  the  Wilke  formula  is  used  (Bird,  Stewart, 
and  Ughtfoot  1960).  This  formulation  is  for  the  low  density  limit,  and  ignores  the  pressure 
dependence  of  the  viscosity.  Fortunately,  for  high  temperatures  the  pressure  dependence  is 
very  small  (Bird,  Stewart,  and  Ughtfoot  1960).  Finally,  a  least-squares  fit  of  the  viscosity  is 
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made  between  300  K  and  3,000  K,  covering  the  temperature  range  of  interest.  The  values  for 
LGP1845  and  LGP1846  are  almost  the  same.  The  expression  used  is 


\i  -  9.53  10-7  7a#“  . 

3.2  Grid.  Figure  4  shows  the  grid  used  in  the  chamber.  The  lower  boundary  is  the  center 
line  of  the  chamber,  and  the  upper  boundary  is  the  wall.  The  solid  lines  represent  the  vector 
grid.  The  grid  is  labeled  from  1  to  nx  in  the  axial  direction  and  from  1  to  nr  in  the  radial 
direction.  The  grid  points  are  evenly  spaced.  The  coordinates  of  the  point  labeled  are 
denoted  by  [xv(i),yv(jj\.  The  axial  velocity  vjij)  and  the  radial  velocity  vy(i,j)  are  defined  at  the 
point  [xv(i),yv(j)\. 

The  vector  grid  splits  the  chamber  into  annular  scalar  control  volumes.  The  scalar  grid  is 
defined  by  the  dotted  lines.  The  scalar  grid  is  labeled  from  1  to  nx-1  in  the  axial  direction  and 
from  1  to  nr-1  in  the  radial  direction.  The  scalar  grid  coordinates  are  defined  by 

xs(/>  [xv(i)  +  xv{i+\)\!2  ys(j)=  [yv(j)  +  yv(/+1)]/2  . 

The  scalar  quantities  are  assumed  to  be  constant  within  each  scalar  volume  (between 
solid  lines).  There  is  a  discontinuity  in  the  scalar  quantities  at  each  solid  line.  Each  scalar 
volume  has  a  volume  (V),  divided  up  into  a  liquid  volume  (VL)  and  a  gas  volume  (VG). 
Similarly,  each  volume  has  a  mass  (M),  divided  into  a  liquid  mass  {ML)  and  a  gas  mass 
The  gas  volume  fraction  e  »  VG/V.  The  gas  has  a  temperature  (T).  The  liquid  is  assumed  to 
be  isothermal,  and  heat  transfer  to  the  liquid  is  ignored.  The  gas  has  density  pc,  and  the 
liquid  density  is  pt.  The  liquid  is  in  the  form  of  drops  with  a  Sauter  mean  diameter  ( d )  and  a 
total  surface  area  (S).  The  pressure  (p)  is  assumed  to  be  the  same  for  both  the  gas  and  the 
liquid. 

The  dotted  lines  split  the  chamber  into  vector  control  volumes.  The  velocities  are 
assumed  to  be  the  same  within  each  vector  control  volume,  and  there  is  a  jump  in  velocity  at 
the  dotted  lines.  Note  that  each  vector  control  volume  is  made  up  of  parts  of  four  scalar 
control  volumes. 
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Figure  4.  Combustion  Chamber  Grid. 


There  are  additional  scalar  control  volumes  around  the  physical  chamber  to  assist  in 
setting  the  boundary  conditions  (see  below). 

3.3  Governing  Equations  -  Combustion.  The  code  takes  a  time  step  in  two  stages.  First, 
a  lagrangian  calculation  is  made.  That  is,  the  grid  is  assumed  to  move  with  the  fluid.  The 
advection  terms  are  considered  in  the  second  stage,  and  the  grid  is  returned  to  the  desired 
location. 

Consider  a  scalar  control  volume  ( V ).  Since  ail  four  vertices  can  have  different  velocities, 
the  Lagrangian  volume  ( V ')  is  no  longer  defined  by  a  simple  rectangle  rotated  about  the 
center  line.  A  general  formula  for  a  quadrilateral  of  rotation  is  required.  Consider  first  one  line 
segment  from  (x„y,)  to  (Xj,y2).  The  volume  under  this  line  segment  after  it  is  rotated  about  the 
center  line  is  given  by 


{it/3)(x2-x,)(y?  +  y,y2+yZ)  . 
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To  find  the  volume  of  the  entire  figure,  apply  this  formula  to  the  four  line  segments  going  in 
the  clockwise  direction.  Bottom  boundaries  will  then  generate  a  negative  volume,  and  the 
total  signed  volume  will  be  correct. 

For  the  Lagrangian  portion  of  the  integration,  there  is  no  motion  across  the  boundaries. 
Each  scalar  control  volumes  can  then  be  integrated  using  the  ordinary  differential  equations 
derived  for  lumped  parameter  models  (Gough  1983;  Coffee  1985, 1988),  with  the  added 
simplification  that  there  are  no  inflow  or  outflow  terms. 

In  the  scalar  control  volumes,  the  densities  are  given  by 

pL «  Ml  Nl  and  pG  =  Mg/Vg. 

The  density  of  the  mixture  is  given  by 


p  »  M/V  -(1  -  to)  pL  +  ta  p a , 


where  eG  -  VG/V.  The  Noble- Abel  equation  of  state  is  used  for  the  gas 


p~paR.T/(1  -b  p6), 


where  R,  is  the  specific  gas  constant  and  b  is  the  covolume.  The  liquid  equation  of  state  is 

p  *<*,/*,>  [(pt/p.>”  -i]. 

where  K ,  is  the  adiabatic  bulk  modulus  at  zero  pressure,  K3  is  the  derivative  of  the  bulk 
modulus  with  respect  to  pressure,  and  pe  is  the  density  at  zero  pressure.  The  pressure  is 
assumed  to  be  the  same  in  both  the  liquid  and  the  gas.  The  speed  of  sound  in  the  liquid  is 
given  by 
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cL  -  ♦  k2p)/ql  , 

where  ge  is  a  conversion  constant,  and  in  the  gas  by 

Co-vkoYp/[Pa(1  -&Pe)J  • 

The  speed  of  sound  in  the  mixture  (assumed  to  be  homogeneous)  is 

c-/(1/p){1/[efi/(pacGa)  +  (1  -e6)/(Ptct2)]}  . 

The  enthaipy  of  the  liquid  is 

/rt-et+p/pt, 

where  the  constant  eL  is  the  chemical  energy  of  the  liquid,  and  the  enthaipy  of  the  gas  is 

ha-cpT  +  bp, 

where  cp  is  the  specific  heat  at  constant  pressure. 

Now  the  governing  ordinary  differential  equations  can  be  written.  Let  m  be  the  combustion 
rate  (g/s).  Then 

dML  /dt » -rh 
and 

dMQ/dt  =  iti  . 

The  total  mass  M  does  not  change.  The  rate  of  change  of  the  volume  is  given  by 

dV/dt  *  (V  -  V)/dt, 
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where  dt  is  the  time  step.  The  rate  of  change  of  the  pressure  is 


dp/dt  «  (p  c2lg0  v)[-dV/dt-m/pL  +  rii/pQ 

*  rhg0{hL  -  h6){y  -  l)/[p6cj(1  -  bpG)]}. 

The  liquid  density  is  found  from  the  equation  of  state.  The  mixture  density  p  =  M/V' .  The 
new  liquid  volume  is  VL *  ML/pL,  and  the  new  gas  volume  is  Va=  V'  -  VL.  The  new  gas 
density  is  then  pG  *  Ma/Va.  The  gas  temperature  is  found  from  the  equation  of  state. 

The  pressure  equation  is  the  only  equation  that  is  integrated  implicitly.  Because  of  the 
large  velocities  in  the  system,  the  volume  can  change  noticeably  in  a  single  time  step.  The 
explicit  integration  may  not  be  accurate  enough.  To  advance  from  time  tn  to  time  tn., ,  the 
derivative  of  pressure  is  found  using  the  known  quantities  at  time  t„.  Then  a  first  estimate  of 
the  pressure  at  time  f„„  is  from  the  explicit  equation 


Pn*i  mP«  +  (dp/dt)„  dt. 


Then  new  values  are  found  for  the  densities  and  sound  speeds.  The  new  Lagrangian  volume 
is  used.  The  same  values  are  used  for  the  rate  of  volume  change,  the  enthalpies,  and  the 
combustion  rate.  Then  the  time  derivative  of  pressure  is  found  at  the  new  time  and 

Pn. I  -  Pn  +  0.5  [  ( dp/dt)n  +  (dp/dt)'.,  ]  dt. 

This  will  converge  after  a  couple  of  iterations. 

To  close  the  system,  the  rate  iti  at  which  the  liquid  droplets  are  combusting  must  be 
known.  The  rate  of  surface  regression  is  assumed  to  be  of  the  form  ApB.  The  rate  of 
combustion  is 

fn  »  p  LS  A  pa  . 
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In  general,  there  will  be  different  sized  droplets  in  each  control  volume.  In  this  code,  only  the 
Sauter  mean  diameter  is  tracked.  If  ail  the  drops  are  assumed  to  have  the  same  diameter 
( d ),  this  is  the  diameter  that  will  preserve  the  surface  area.  For  one  droplet, 

V,mK  d*/6  and  S,  *  jc  d*. 

So  the  total  liquid  surface  is  given  by 

S*  6  VL/d, 

and  the  total  number  of  drops  is 

N  m  Vl/V1  *  VL  6/it  d3 . 

The  quantity  of  interest  is  the  surface  area,  S.  This  can  also  be  written  as 

S-A/icd*. 

Then 

dS/dt  -  N  k  2  d  dd/dt 

-  Nk  2 d{  -  2 A  pa) 

-  -  (Vt  6/%  d3)  K2d{2fnlpLS) 

-- VL(24/d2)(fn/pJ(d/6VL ) 

■  -  4  fti  /  dpL  . 

The  new  Sauter  mean  diameter  can  then  be  found  from  the  surface  area 

If  there  is  no  influx  or  outflux,  mass  and  energy  conservation  should  hold.  The  energy  in  a 
scalar  control  volume  is  given  by 

I  *  ML  eL  +  MG  cvT  +  M  v*  /  2  g0 , 
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where  v  is  the  average  magnitude  of  the  velocity  in  the  control  volume.  The  last  term  is,  in 
practice,  found  in  four  parts.  The  magnitude  of  the  velocity  at  each  vertex  of  the  control 
volume  is  found,  and  this  is  multiplied  by  the  appropriate  mass. 

3.4  Time  Step.  Now  we  can  consider  the  procedure  for  advancing  a  time  step.  Consider 
the  axial  velocity  equation.  The  integration  is  performed  over  a  vector  control  volume 
(bounded  by  dotted  lines).  The  velocity  vji.j)  is  assumed  constant  over  the  control  volume, 
so  it  can  be  taken  outside  the  integral.  The  density  can  take  on  four  different  values,  since 
the  vector  control  volume  is  made  up  of  parts  of  four  scalar  control  volumes.  Let  Vt(J)  be  half 
the  volume  of  the  lower  section  of  the  scalar  volume  Vfj)  and  VJJ)  be  half  the  volume  of  the 
upper  half.  Then,  if  Mv  is  the  mass  in  the  vector  control  volume, 

Mv  -  Jyp  dV  -  p(/.y)  *  V,  (j)  ♦  p(/- 1 ,/)  *  V,  (/) 

♦  p(U-i)*v„(y-i)  +  p(/-i,/-i  )*vu{j-t ) . 

The  advection  terms  (first  two  terms  on  the  right)  disappear  for  the  Lagrangian  step.  The 
momemum  equations  have  actually  been  given  for  a  fixed  grid,  in  general,  the  advection 
terms  are  given  by  the  momentum  p  v,  times  the  velocity  across  the  control  volume  boundary. 
For  a  Eulerian  calculation  (fixed  grid),  the  velocity  across  the  axial  boundaries  is  vx ,  and  the 
velocity  across  the  radial  boundaries  is  v, .  For  a  Lagrangian  calculation,  there  is  no  flow 
across  the  boundary. 

There  are  pressure  terms  on  the  right  and  left  boundaries.  Let  >4,(/)  be  the  lower  axial 
area  of  a  scalar  control  volume  and  Au(j) )  be  the  upper  area.  Then 

d(Mvv,)( i,j)/d f  -  p( / - 1 ,/) *A, ( j )  +  p(i - 1  ,j - 1 ) mAu{j - 1 ) 

-p{/./)*A(/)  -p(/./-1).^(y-l)  +  ... 
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The  viscosity  term  is  more  complicated.  Let 


dvl „(/,/)  «  2vjdx  on  left  boundary  of  the  vector  volume  (/,/) 

dv\n(i,j)  •  dvy/dy  on  left  boundary  of  the  vector  volume  (/,/) 

*  0.5[vy(/,/+1)  -  vy(i,j- 1)  +  vy(/-1,/+1) 

-  v,(/-1J-1)]/[yv</*1>  -yv(y-l)]; 

dvbv(l,fl  ■  dV'/tjx  on  bottom  boundary  of  the  vector  volume  {/.;) 

-  [vAl.j)  -  M/J-l  )]/[yv(y)  -  yv{J- 1 )]; 

dvb^fij)  »  dvyldx  on  bottom  boundary  of  the  vector  volume  (Ay) 

-  0.5*[vy(/+1,y)  -  v^i-fA)  +  vy(/+1,y-1) 

-  *V(/-1.y-1)]/[xv(/+1)  -  xv {I  —  1 ) ] . 

The  axial  viscosity  term  becomes 

3(  Mvvt)  (i,j)/dt  -  ...  -[4/3  dvl„{i,j )  -  212  dvlyy(i,j)] 

[n(/'-i.y)A(y)  ♦  |i(/-1,y-1  Mu(y-1 )] 

+  [4/3 dvlxx{i+ 1  ,y)  -  2/3dW,y(/+1  ,y)] 

[M'.y)A(y)  ♦  h(/.7-i)A,(7-i)]- 
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There  is  an  added  complication  in  the  radial  viscosity  term.  On  each  boundary,  we  want  the 
viscosity  term  times  the  radial  unit  vector  times  the  normal  vector.  But  since  the  area  of  a 
control  volume  increases  in  the  radial  direction,  this  term  is  not  zero  on  the  axiai  sides  of  the 
control  volume.  Under  the  assumption  that  the  viscosity  term  varies  linearly  from  the  bottom 
to  the  top  boundary,  a  little  algebra  shows  that  the  total  contribution  is  the  difference 
between  the  two  sides  times  the  area  halfway  from  the  bottom  to  the  top  side.  Call  this  area 
Am(j)  (area  in  the  middle  of  the  vector  control  volume  j  or  at  the  bottom  of  the  scalar  control 
volume  J).  Then 


B{Mvvg){i,j)/dt  =  ...  -  I dvbxy(i,j)  +  dvbyx{i,j)]  0.5 

ln(/-1.y-1)  +H(/,/-1)]4n(/) 

♦  +  cfvfrr,(/,y+i)]o.5 

M/-1.7)  +\L(i.j)]AJj^)  . 

For  the  radial  momentum  equation,  we  also  need  to  define 

dvlxy(iJ)  *  dvx/dy  on  left  boundary  of  the  vector  volume  (/,;') 
.0.5[v„(/,/.1>  -  ♦  v,(/-1,/.1) 

-  *i(/- 1 1  )]'[yv(J* i )  -  yvu- 1 )] ; 

*  dVy/dx  on  left  boundary  of  the  vector  volume  (/./) 

-  [Vy(ij)  -  vyV-l.j)]/[xv{i)  -xv(i- 1)]; 

dvb^i.j)  «  dvy/dy  on  bottom  boundary  of  the  vector  volume  (/,/) 

■  [ vy(i’j)  -  vy(i.j-l)]l[yv{j)  -  yv(j- 1 )]; 
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dvbu(i,j)  ■  dvM/dx  on  bottom  boundary  of  the  vector  volume  (/,/) 

-  0.5*[vx(/+1,y)  -  v,</-iy)  +  vx(i+ 1  f] ~  1) 
-^(/-i,y-i)]/[xv(/+i)-xv(/-i)]. 

a  ( Mvvy)  ( /,y)/3f  -  o.5[p( /- 1  ,y- 1 )  +  p(  /,y- 1) 

-p(/-i.yj  -p(/./)]^m(y-i) 

-  (w) +  <H,(/y)] 

[n(/-i.y)A(7)  +  n(/-i,y-i)>uy-i )] 

♦  I<Hy(/+1J) 

[m.(/,/)  a,  (7)  +ii(/.y-i)^(y-i)] 

-  [4/3  dvbyy(i,j)  -  2/3  cfvb„(  ij) ]  0.5 

IU(/-1.7-1)  ♦H(/J-1)1A.(/) 

+  [4/3cW>yr(/,y+1 )  -  2/3  dvbxx(i,j+ 1 )  ]  0.5 

[n(/-i.y)  +  n(/,y)]^m(y+D. 

TThe  new  value  of  mass  times  velocity  is  found  by  taking  the  old  value  plus  the  time  derivative 
times  the  time  step.  The  quantities  Mvvx  and  Mvvy  are  kept  in  vectors.  The  actual  velocities  vx 
and  vy  are  also  updated  at  this  time. 

Next  consider  the  coordinates  [xv{i),  yv(j)]  of  a  vector  grid  point  (intersection  of  solid 
lines).  Let  xvl(itj)  and  yvl(i,j)  be  the  corresponding  coordinates  at  the  end  of  the  Lagrangian 
part  of  the  time  step  (this  grid  is  not  orthogonal).  Then 

xv/(/,y)  -  xv(/)  +  vx(i,j)dt ; 
yvl{i,j)  •  yv(j)  *  vy(i,j)dt. 
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The  Lagrangian  volume  V1  of  the  scalar  control  volumes  can  now  be  found  (see  Equation  5). 
Then  the  equations  in  Section  B  can  be  applied,  and  ail  the  scalar  quantities  can  be  updated. 

It  might  be  more  consistent  to  use  the  velocities  at  the  end  of  the  last  time  step,  rather 
than  the  velocities  calculated  after  the  pressure  step.  However,  the  integration  is  more  stable 
using  the  updated  velocities.  The  flow  is  kept  current  with  the  present  pressures. 

Next,  the  grid  is  updated.  In  general,  the  left  wall  of  the  combustion  chamber  is  allowed  to 
move  at  an  arbitrary  speed,  vleft  Then  the  grid  velocity  at  the  left  vg(1)  =  vleft.  The  right  wall 
is  kept  fixed  [vg(nx)-  0.0).  The  axial  grid  is  kept  evenly  spaced,  so  the  grid  velocity  is 
assumed  to  be  a  linear  function  of  position.  The  grid  velocities  at  the  other  values  are  found 
by  interpolation.  Then  the  new  position  of  the  axial  vector  grid  is 

xv(i)  ■  xv(/)  +  vg(l)dt . 

The  radial  grid  positions  are  unchanged. 

Now  the  advection  terms  for  the  scalar  quantities  can  be  found.  Conceptually,  the 
boundaries  of  the  Lagrangian  control  volumes  are  moved  to  the  boundaries  of  the  new  control 
volumes.  The  amount  of  material  the  boundary  passes  through  is  convected  from  one  volume 
to  the  next.  For  instance,  consider  the  left-hand  boundary  of  the  scalar  control  volume  (Ay). 
Assume  that  both  Lagrangian  points  are  to  the  right  of  the  final  boundary  points.  Compute  the 
volume  of  the  rectangle  bounded  by  the  points  [xv(/),  yv{j+1)],  xvf(i,j+1),  yvl{i,j+1),  [xvl(ij), 
yvl{ij)],  and  [xv{i),  yv(J) ],  integrating  in  this  order.  For  the  assumed  case,  this  is  a  positive 
number,  which  implies  that  the  flow  is  from  the  previous  control  volume  (M,y)  to  the  present 
volume  (/,/).  Upwind  differencing  is  used,  so  the  value  of  any  quantity  on  the  boundary  is  set 
equal  to  the  value  in  the  volume  (/-/,/).  In  this  case,  the  upwind  differencing  is  physically 
correct,  since  quantities  will  be  moved  from  the  volume  (i-1,j)  to  the  volume  (Ay).  If,  instead, 
both  Lagrangian  points  are  to  the  left  of  the  final  grid  boundary,  the  volume  is  negative. 
Material  will  be  convected  from  the  volume  (Ay)  to  the  volume  (M.y).  Again,  using  upwind 
differencing,  the  correct  values  will  be  used.  The  ambiguous  case  is  when  one  Lagrangian 
point  is  to  the  left  of  the  boundary  and  the  other  point  is  to  the  right.  Some  material  will  be 
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converted  to  the  right  and  some  material  to  the  left.  The  computation  above  will  give  the 
proper  signed  volume  change  for  the  volume  (/,/).  However,  upwind  differencing  is  no  longer 
exactly  correct.  The  upwind  differencing  is  used  because  it  is  stable  even  in  the  cases  when 
it  is  not  strictly  correct,  and  for  this  formulation  it  is  more  often  correct  than  other  simple 
differencing  schemes. 

Consider  first  the  gas  mass.  Let  VI  be  the  volume  between  the  left  boundary  of  the 
Lagrangian  volume  and  the  left  boundary  of  the  final  boundary.  Then  set 


fluxl  »  [ /WG(/- 1 ,/) /  V' (  / - 1 ,/) ]  Vl(i,j) ,  if  W(/./)>0 
[Ma(i,j)IV'(i.j)]vi(i,j),  if  VI(IJ)< 0  . 


The  flux  through  the  other  three  boundaries  is  computed  the  same  way.  Then  the  new  gas 
mass  is 

Ma(i,j)  •  Ma(i,j)  +  ( fluxl  +  fluxr  +  fluxb  +  fluxt) . 

For  the  liquid  mass,  the  analogous  equations  are 

fluxl  -  [  Ml(  I  - 1 ,/)  /  V  ( /  - 1 ,» ]  VI ( /,/) ,  if  W(  /,/)  >  0 
[ML(i,j)IV'(i,j)]vi(i,j),  if  W(/,/)<0 


and 


ML{i,j)  ■  ML(i,j)  +  (fluxl  +  f/uxr  +  fluxb  +  fluxt). 
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For  the  energy  equation 


fluxl  -  [ Me(i- 1 ./)/  V'(i- 1  ,j)]  7(/- 1  ,i)  Vl(i,j) ,  if  W(/./)  >  0 

[wt(/.y)/ W(/.y)]r(/.y)  w(/,y).  if  w(/.y)<o 


and 


( M0T)(i,j )  -  (MaT)(i,j)  +  { fluxl  +  fluxr  +  fluxb  +  tfuxf). 


There  is  also  an  equation  for  the  surface  area. 


fluxl  -  [S(/- 1  ,y>/ w (/- 1  .y) ]T(i- 1 ,/)  w(/,y) ,  if  w(/,yj > o 
ls(ij)/v'(i,j)]T(i,j)  vi(ij),  if  w(/,yj<o 

and 

S(/,y)  •  S(/,y)  +  ( fluxl  +  fluxr  +  fluxb  +  fluxt) . 

A  similar  procedure  is  used  for  the  vector  quantities.  The  vector  control  volumes  are 
bounded  by  the  scalar  points.  The  boundary  locations  of  the  vector  control  volumes  are  given 
by 

xsi(i.j)  - [xv/(/,y+i )  +  xw(/+ i,7+i )  +  xv/(/+i,y)  +  xv/(/,y)]/4 
ys/(/,y)  -  [yv/(/,y+1 )  +  yW(/+ 1,7+1)  +  yv/(/+1  J)  +  yv/( /.y) ]/4  . 


Let  V/(/,y)  now  represent  the  volume  change  at  the  left  boundary  of  the  vector  control  volume 
from  the  Lagrange  to  the  final  position.  The  density  on  the  left  boundary  is  approximated  by 
an  area  weighted  average 
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p/(/.y)  -[p(/-i.y-i)40-i.y-i)  +p(/-i./)a(/-i./)]/ 

[>A„(/-i.y-i )  ♦  a(/-i./)]. 

For  the  axial  velocity, 

fluxl  -  p /{/,/)  v,(/-1 ,J)  W,  W(/,y)  >  0 

p/(/./K(/./)  W.  Vl(i,j)<0 

and 

(Myvx)(i,j)  •  (Mvvx)(i,j)  +  ( fluxl  +  fluxr  +  fluxb  +  fluxt). 

Similarly,  the  radial  velocity  is  given  by 

fluxl  -  p/(/,7)  v,(/-1,y)  W  W(/,y)>0 

pl{IJ)vy(i.j)VI  Vl(i,j)<0 

and 

(Mvvy)(i,j)  ■  ( Mvvy){i,j )  +  ( fluxl  +  fluxr  +  fluxb  +  fluxt) . 

The  formulas  need  to  be  modified  for  the  vector  control  volumes  along  the  centerline.  In 
this  case,  the  standard  form  of  the  vector  control  volume  would  extend  below  the  centerline. 
However,  since  a  volume  of  rotation  is  being  considered,  this  does  not  add  anything  to  the 
vector  control  volume.  So  for  a  velocity  at  the  centerline,  the  vector  control  volume  is 
considered  to  be  only  the  part  above  the  centerline.  So  only  the  part  actually  in  the  physical 
chamber  is  considered.  Special  logic  is  then  required  for  the  change  in  volume.  For 
consistency,  the  vector  control  volumes  along  the  top  boundary  are  treated  the  same  way. 
That  is,  the  vector  control  volume  does  not  extend  past  the  physical  boundary. 

Most  of  the  other  quantities  can  now  be  calculated  in  a  straightforward  manner.  The 
exception  is  the  pressure  term.  By  assumption,  the  pressure  in  the  new  control  volume  ( /,/ )  is 
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uniform  and  equal  for  both  the  gas  and  the  liquid.  The  new  pressure  is  somewhere  between 
the  pressure  in  the  Lagrangian  volume  and  the  pressure(s)  of  the  fluid  convected  into  the  final 
scalar  volume.  But  as  the  pressure  changes,  both  the  liquid  and  gas  volumes  change.  The 
constraint  is 


PL-Po  =  0.0 

From  the  gas  and  liquid  equations  of  state 

(K1//C,)[(pt/p.)«  -  l]  -  psfl.77(1  -  Op.)  •  0 

Rearranging,  a  function  of  the  liquid  volume  can  be  defined 

f(VJ  -  {K,IK,)(V-  VL  -  bUc)[{ML/p.r  -  Vi."] 
-  R,MaT  VL*2  .  0  . 


The  equation  is  transcendental  and  is  solved  by  the  Newton-Raphson  iteration.  The  derivative 

-  -(K,/K,)[(Mt/p„)”  -  Vi”] 

-K,(V-VL-bMa)  V?" 

-  R,MaTKtV”-' 


and 

VL.VL-f(VL)/f'(VL). 

The  time  step  is  determined  by  the  Courant-Friedrichs-Lewy  condition  (McBratney  1980),  as 
modified  for  compressible  flow 
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dt£dx/[c  +  vx] 
dt£dy/[c  +  vy\  . 


That  is,  a  signal  cannot  propagate  farther  than  a  cell  length  during  one  time  step.  The  code 
loops  over  all  the  scalar  control  volumes  and  finds  the  minimum  dt  by  the  above  criteria.  The 
velocity  in  the  cell  is  taken  to  be  the  average  of  the  four  corner  velocities.  The  time  step  is 
then  multiplied  by  a  user-supplied  safety  factor  (usually  0.5).  This  is  further  modified  so  that 
the  output  times  are  hit  exactly. 

3.5  Boundary  Conditions.  An  extra  set  of  scalar  control  volumes  are  placed  around  the 
physical  chamber.  These  volumes  are  the  same  size  as  the  scalar  control  volumes  just  inside 
the  physical  boundary,  except  at  the  top  boundary. 

Consider  the  case  where  the  chamber  is  completely  closed  (no  inflow  or  outflow).  The 
physically  correct  boundary  condition  on  the  walls  is  a  no  slip  condition.  That  is,  the  velocities 
on  the  boundary  are  zero.  This  is  included  as  an  option  in  the  code. 

For  this  case,  consider  a  scalar  control  volume  next  to  the  physical  boundary.  To 
integrate  the  scalar  quantities,  the  four  sets  of  velocities  at  the  comem  of  the  control  volume 
are  required.  These  are  known.  Values  for  the  four  scalar  control  volumes  adjacent  to  the 
volume  of  interest  are  also  required.  This  includes  one  scalar  control  volume  outside  the 
physical  chamber.  However,  there  is  no  convection  across  the  boundary.  So  while  these 
control  volumes  are  put  in  to  avoid  having  a  special  code  for  the  boundaries,  it  is  irrelevant 
what  values  are  assigned  to  the  extra  volumes.  Consider  the  vector  control  volume  nearest 
the  boundary.  The  values  in  the  four  control  volumes  that  make  up  the  vector  control  volume 
as  well  as  the  eight  velocities  in  a  box  around  the  point  of  interest  are  required.  These  are  all 
available. 

At  the  centerline,  the  radial  velocities  must  be  zero,  so  there  is  no  mass  flow  across  the 
centerline.  The  axial  velocity  under  the  boundary  (j  =  0)  is  given  the  same  value  as  the  axial 
velocity  over  the  boundary  (j  =  1).  The  radial  velocity  is  reversed  (flow  up  becomes  flow 
down). 


25 


For  the  high  pressures  and  velocities  in  a  gun  simulation,  the  boundary  layers  will  be  very 
small,  and  it  is  impractical  to  resolve  these  in  detail.  Instead,  a  slip  condition  is  normally  used. 
It  is  assumed  that  the  boundary  layer  is  negligibly  small  compared  to  the  grid  spacing,  and  the 
boundary  layer  is  ignored.  Instead,  the  tangential  velocity  is  calculated  using  the  standard  set 
of  equations.  The  normal  velocity  is,  of  course,  still  zero.  This  has  no  effect  on  the  scalar 
quantities,  since  the  flux  through  the  boundary  is  still  zero.  However,  the  tangential  velocities 
at  the  boundaries  must  be  calculated. 

Consider  a  point  on  the  top  boundary.  The  scalar  quantities  just  outside  the  boundary  are 
given  the  same  values  as  the  scalars  just  inside  the  boundary.  However,  the  axial  areas  are 
set  to  zero.  Then  the  pressure  equation  can  be  used  unchanged.  To  obtain  the  proper 
velocity  derivative  terms,  extrapolation  conditions  are  used  for  the  velocities  —that  is, 

vx(l,nr+ 1 )  »  vx(i,nr)  +  [vx(i,nr)  -  vx(i,nr- 1 )] 
vy(i,nr+ 1 )  -  vy(i,nr)  +  [ vy(i,nr )  -  vy{i,nr- 1 )]  . 


On  the  right  and  left  boundaries  analogous  conditions  are  used. 

The  code  allows  for  injection  of  liquid  through  an  annulus  in  the  left  boundary.  The  upper 
and  lower  boundaries  yu  and  y,  are  input,  as  well  as  the  injection  rate  m*,  in  g/s,  and  the 
injected  droplet  diameter  d The  injection  area  is  denoted  by  A^.  Consider  the  scalar  control 
volumes  which  are  totally  or  partially  inside  the  injection  annulus.  The  corresponding  scalar 
control  volumes  just  to  the  left  of  the  boundary  are  assumed  to  be  pure  liquid  with  the  density 
of  the  liquid  just  inside  the  boundary  but  made  of  droplets  with  the  diameter  d f*, .  Pure  liquid 
cannot  really  be  made  up  of  droplets,  but  this  is  a  way  of  simulating  the  assumption  that  the 
liquid  breaks  up  into  droplets  of  the  specified  size  as  soon  as  it  enters  the  chamber.  Now 
consider  the  vector  points  between  yu  and  yr  Let  rJJ)  be  the  ratio  of  the  area  of  the  left 
boundary  of  the  vector  control  volume  that  is  in  the  injection  annulus  to  the  entire  control 
volume  area.  That  is,  the  vector  control  volume  might  only  be  partially  in  the  annulus.  Then 

v,(1.y)  «  vleft  +  rm(j) ( rhlnl Aln(j) ) /  0.5  [pt(1 ./)  +  pjl.y'-l )]  . 


26 


The  other  boundary  quantities  are  found  as  before. 

There  is  a  problem  if  the  injection  area  is  smail  compared  to  the  ghd.  The  amount  of 
liquid  injected  depends  on  the  Lagrange  volume  computed.  The  Lagrange  volume  calculation 
assumes  that  the  gas  velocity  varies  linearly  from  one  vector  grid  point  to  the  next.  For  the 
injection  process,  however,  there  is  a  physical  boundary  where  the  velocity  should  suddenly 
change  to  zero.  To  obtain  the  proper  injection  rate  even  for  a  coarse  grid,  special  logic  is 
used  to  compute  the  volume  change  over  the  injection  area.  Tht>  volume  change  at  the  left 
boundary  of  a  scalar  control  volume  is  the  sum  of  the  axial  velocities  at  the  boundary  times 
the  appropriate  vector  areas  times  the  time  step.  The  injection  area  is  allowed  to  vary  with 
time. 

There  is  also  a  primer  injection  model.  In  the  gun,  the  primer  pressurizes  the  combustion 
chamber  to  start  the  firing  cycle.  The  actual  primer  is  injected  through  a  circular  opening  and 
can  only  be  modeled  by  a  three-dimensional  code.  However,  the  details  of  the  primer 
injection  are  not  of  interest,  and  the  requirement  Is  just  to  bring  the  chamber  up  to  about  the 
correct  pressure.  The  model  used  in  the  lumped  parameter  code  is  extended  to  the  present 
case.  The  primer  is  assumed  to  inject  a  specified  mass  of  hot  gas  Mp  uniformly  over  a  time 
period  tp.  The  injection  rate  is  then  rfy  «  Mp/tp.  A  top  and  bottom  boundary  in  the  right  wall  is 
specified,  and  the  injection  area  Ap  is  calculated.  Normally  about  50%  of  the  primer  energy  is 
lost  during  the  injection  process.  So  an  energy  loss  term  fp  is  also  input. 

* 

During  the  primer  injection,  the  scalar  control  volumes  outside  of  the  injection  area  are 
assumed  to  be  pure  gas.  The  density  is  more  or  less  arbitrarily  set  to  0.2.  The  primer  is 
assumed  to  be  the  same  material  as  the  liquid  propellant  The  temperature  of  the  gas  is  the 
energy  term  fp  times  the  chemical  energy  of  the  propellant  divided  by  the  specified  heat  at 
constant  volume.  Consider  all  the  vector  points  in  the  annulus.  Let  rp  be  the  ratio  of  the  area 
of  the  right  boundary  of  the  vector  control  volume  that  is  in  the  injection  annulus  to  the  entire 
control  volume  area.  Then 

^(nx.j)  -  -  rp(j)[mp/Ap(j)]/ 0.5  [pG(nx,j)  +  p0(nx,j- 1 )]  . 

After  the  primer  has  been  injected,  this  boundary  is  treated  like  a  standard  wall. 
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Finally,  there  is  the  flow  into  the  gun  tube.  This  is  a  circular  opening  at  the  centerline  of 
the  chamber.  The  gun  tube  must  not  overlap  the  primer  injection  annulus.  The  radial  grid  is 
adjusted  so  that  the  top  of  the  gun  tube  always  coincides  with  a  vector  grid  point.  Since  the 
formulation  of  the  equations  does  not  depend  on  the  radial  grid  being  evenly  spaced,  this 
does  not  cause  any  problems.  The  gun  tube  can  be  modeled  in  detail  (see  below),  or  this 
can  be  viewed  just  as  an  outflow  boundary. 

If  the  detailed  gun  tube  model  is  not  used,  the  pressure  p,  just  outside  the  chamber  must 
be  specified.  This  is  the  pressure  in  the  scalar  control  volumes  just  outside  the  chamber  in 
the  tube  region.  The  flow  through  the  boundary  is  assumed  to  be  isentropic.  For  a 
Noble-Abel  equation  of  state,  the  process  equation  is 

p(  1/p  -b)1  -  constant . 

Then  the  gas  density  can  be  found  from  the  process  equation  and  the  gas  temperature  from 
the  equation  of  state.  The  mass  of  liquid  just  outside  the  chamber  is  assumed  to  b“  the  same 
as  the  mass  of  liquid  just  inside  the  chamber.  The  liquid  density  is  found  from  the  liquid 
equation  of  state.  The  velocity  just  outside  the  chamber  is  found  by  extrapolation. 

There  is  a  sharp  comer  at  the  top  of  the  gun  tube  region.  The  top  right  comer  of  the 
vector  control  volume  for  this  point  is  not  actually  in  the  fluid.  There  is  flow  out  of  the 
chamber  through  the  bottom  part  of  the  vector  volume,  but  the  top  part  of  tl}e  vector  volume 
still  includes  the  chamber  wall.  So  the  axial  velocity  at  this  comer  is  considered  only  to  apply 
to  the  bottom  part  of  the  vector  volume.  Similarly,  the  radial  velocity  only  applies  to  the 
left-hand  part  of  the  vector  volume.  Special  logic  is  put  into  the  integrator  to  handle  this  one 
special  case. 

For  some  problems,  there  may  be  flow  into  the  chamber.  The  above  boundary  conditions 
are  not  consistent  for  inflow,  since  extrapolation  should  be  done  in  the  upwind  direction. 
Generally,  the  conditions  outside  the  chamber  are  not  known  in  detail.  So  to  preserve 
stability,  back  flow  into  the  chamber  is  not  allowed.  If  the  flow  reverses,  the  exit  velocity  is 
arbitrarily  set  to  zero. 
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For  guns,  the  flow  out  of  the  chamber  will,  in  general,  be  subsonic.  There  is  a  test  fixture 
at  Sandia  Laboratories  (see  below)  where  the  chamber  empties  into  atmospheric  pressure  and 
the  flow  is  choked.  A  burst  disc  prevents  flow  until  the  chamber  is  properly  pressurized. 

To  model  this  case,  the  outflow  is  set  to  zero  until  the  pressure  at  the  right  centerline  goes 
over  the  burst  disc  pressure.  Then  the  assumption  is  made  that  choked  flow  is 
instantaneously  established. 

For  an  ideal  gas,  there  are  algebraic  formulas  for  choked  flow  (Fox  and  McDonald  1985). 
Let  M  be  the  mach  number  of  the  gas  on  the  centerline. 

Mm  ^(nx-l.l  )/c(nx-1,1 )  . 

The  stagnation  pressure,  p^  is  the  pressure  of  the  gas  if  the  velocity  is  isentropically  reduced 
to  zero. 

pa  •Pinx  -  1,1 )  [l  ♦  (y-  1  )Ma  /  2],t/<7_1>1  . 

The  critical  pressure  is  the  pressure  at  which  the  Mach  number  is  unity. 

Pe»Po[(Y+1)/2]t'1r/<T_nl. 

* 

For  choked  flow,  the  pressure  at  the  throat  is  the  critical  pressure,  and  the  velocity  is  the  local 
speed  of  sound.  The  stagnation  enthalpy  of  the  fluid  is  defined  as 

h0  -  h  *  v2/2  , 

and  this  is  constant  throughout  an  adiabatic  flow  field. 


For  the  present  two-phase  mixture,  where  the  gas  follows  the  Noble-Able  equation  of 
state,  an  analytic  solution  is  not  possible.  Instead,  an  iterative  procedure  is  used.  The  initial 
guess  for  the  value  of  pc  is  found  using  the  above  ideal  gas  equations.  The  corresponding 


gas  density  is  found  from  the  process  equation,  and  the  gas  temperature  from  the  equation  of 
state.  The  mass  fraction  *  MJ(Ma  +  MJ  is  assumed  to  be  the  same  as  just  inside  the 
chamber.  The  liquid  density  is  found  from  the  liquid  equation  of  state.  The  enthalpy  of  the 
fluid  in  the  throat  is  given  by 

•  (CpTe  *  t>Pc)  *ug  *  (et+Pc/pu)0  -e**). 

The  throat  velocity  is  found  from  the  stagnation  enthalpy  equation.  The  throat  speed  of  sound 
is  found  from  the  usual  equations.  The  difference  between  the  speed  of  sound  and  the  fluid 
velocity  is  computed.  If  this  is  less  than  10'7,  the  iteration  is  concluded.  Otherwise,  the  critical 
pressure  is  modified  by 

P.-P.I1 

and  the  calculation  is  redone.  Once  the  iteration  concludes,  the  conditions  outside  the 
chamber  are  set  equal  to  the  critical  conditions,  and  the  velocity  at  the  exit  is  the  speed  of 
sound  at  the  critical  conditions.  However,  in  the  code,  the  velocity  acts  on  the  fluid  just  inside 
the  chamber.  So  the  critical  sound  speed  is  multiplied  by  the  ratio  of  the  critical  density  and 
the  density  just  inside  the  chamber.  This  gives  the  appropriate  mass  flow  out  (density  times 
velocity  times  area).  For  simplicity,  the  calculation  is  only  done  for  the  centerline  values,  and 
the  exit  velocity  is  assumed  to  be  uniform.  The  radial  velocities  at  the  exit  are  set  to  zero. 

3.6  Combustion  Rate.  All  the  problems  considered  in  this  report  involve  LGP1846.  The 
surface  regression  rate  has  been  measured  as 

1.64  ff 103  0  <  p  <  60 

by  McBratney  (1980,  1981).  McBratney  observed  some  indication  of  a  slope  break  in  the 
burning  rate  above  60  MPa.  More  recent  work  by  Oberle  and  Wren  (1990)  indicates  a 
regression  rate  of 


0.000577  p*  009  100  <p  <200  . 
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In  the  code,  a  two-part  burning  rate  is  used,  with  the  McBratney  rate  used  under  67  MPa  and 
the  Oberie  and  Wren  rate  above  67  MPa. 

Due  to  the  rapid  burning  rate,  very  high  pressures  can  be  generated  near  the  injector.  To 
simplify  the  numerics,  the  burning  rate  is  assumed  to  reach  a  maximum  at  1 ,000  MPa.  That 
Is,  for  pressures  above  1 ,000  MPa,  the  burning  rate  is  set  equal  to  the  value  at  1 ,000  MPa.  In 
actuality,  the  burning  rate  is  expected  to  level  off  at  some  value,  but  the  point  where  this  will 
happen  is  not  known. 

To  set  up  acoustic  pressure  oscillations,  there  must  be  some  liquid  accumulation  in  the 
chamber.  Initial  combustion  generates  a  pressure  wave.  The  pressure  wave  bounces  off  the 
chamber  wall  and  returns  to  the  combustion  region.  The  higher  pressure  increases  the 
combustion  rate,  which  further  increases  the  pressure.  The  pressure  wave  then  becomes 
steadily  larger  until  nonlinear  effects  come  into  play.  From  studies  of  liquid  propellant  rockets, 
it  is  known  that  acoustic  oscillations  can  only  be  set  up  if  the  pressure  exponent  of  the 
regression  rate  is  large  enough  (generally  greater  than  one)  (Harrje  1972).  In  fact,  if  the 
McBratney  rate  alone  is  used,  the  code  will  not  generate  pressure  waves  under  any 
circumstances.  For  infinitesimal  pressure  waves,  the  natural  frequencies  of  a  chamber  can  be 
solved  for  analytically  (acoustic  modes)  (Harrje  1972).  Even  for  large  pressure  waves,  the 
observed  frequencies  are  usually  close  to  the  acoustic  modes. 

4.  VALIDATION 


Some  simple  test  problems  were  set  up.  The  various  options  were  tested  with  a  closed 
chamber  (no  outflow)  and  then  with  an  open  chamber. 

4.1  Closed  Chamber.  First,  considr  i  closed  chamber  5.0  cm  long  and  5.0  cm  in 
diameter.  The  chamber  was  filled  with  a  uniform  gas  mixture  at  100  MPa  and  2,000  K.  The 
code  was  run  for  several  milliseconds.  The  velocities  become  nonzero  because  of  round-off 
error.  However,  the  velocities  stay  negligibly  small,  and  the  pressures  and  temperatures  do 
not  change  noticeably.  This  indicates  that  the  code  is  stable. 
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Next,  a  10-cm-long  chamber  is  considered.  The  left  side  of  the  chamber  is  filled  with  gas 
at  50  MPa  and  the  right  side  at  100  MPa.  At  time  zero,  a  rarefaction  wave  propagates  to  the 
right  and  a  shock  wave  to  the  left.  This  case  can  be  solved  analytically  (Hughes  and  Brighton 
1967;  Band  1960).  The  rarefaction  wave  spreads  out  as  it  propagates.  The  front  of  the  wave 
moves  at  the  speed  of  sound,  and  the  back  of  the  wave  moves  at  the  sound  speed  minus  the 
gas  velocity.  The  shock  wave,  on  the  other  hand,  will  become  steeper. 

With  the  no  slip  boundary  conditions,  this  reduces  to  a  one-dimensional  problem.  In  the 
code,  all  the  quantities  are,  in  fact,  uniform  in  the  radial  direction.  Figure  5  compares  the 
theoretical  solution  with  four  different  axial  grids  after  a  short  time  period.  Even  with  the 
coarse  grid,  the  solution  is  reasonable.  The  fine  grid  is  very  close  to  the  theoretical  solution. 

There  is  one  problem  here.  The  pressure  between  the  shock  wave  and  the  rarefaction 
wave  should  be  flat.  Instead,  there  is  a  rise  in  pressure  just  to  the  right  of  the  shock  (Gibbs 
phenomenon).  This  is  an  artifact  of  the  numerical  procedure  and  cannot  be  eliminated  by 
grid  refinement  There  are  techniques  for  eliminating  this  phenomenon.  Artificial  viscosity  can 
be  added  to  the  equations.  However,  this  will  smear  out  all  the  peaks  in  the  solution, 
including  the  pressure  waves  of  interest.  A  more  sophisticated  procedure  is  flux-corrected 
transport  (Peyret  and  Taylor  1983;  Boris  and  Book  1976).  This  Is  a  procedure  for  diffusing  out 
numerically  introduced  peaks  while  having  minimal  effect  on  the  rest  of  the  solution. 
Unfortunately,  for  the  problems  of  interest,  some  high  frequency  oscillations  will  only  be  a  few 
grid  points  in  size,  and  the  flux-corrected  transport  will  also  tend  to  damp  these  physical 
peaks.  Flux-corrected  transport  has  been  implemented  as  an  option  in  the  code  but  will 
usually  not  be  used.  For  a  fine  enough  grid,  essentially  the  same  solution  is  obtained  with 
and  without  flux-corrected  transport. 

Next,  a  uniform  mixture  of  gas  and  liquid  is  put  in  the  chamber.  The  liquid  is  assumed  to 
be  LGP1846,  which  has  a  flame  temperature  of  (eL/cJ  =  2,468.5  K.  A  total  of  5  g  of  liquid  is 
distributed  in  the  chamber,  with  a  droplet  diameter  of  100  pm.  The  initial  gas  temperature  is 
set  to  the  flame  temperature.  At  the  end  of  the  combustion,  there  is  a  uniform  gas  mixture  at 
the  flame  temperature.  So  the  combustion  procedure  works  properly.  Both  mass  and  energy 
are  conserved. 
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The  liquid  was  then  axially  put  into  the  first  centimeter  of  distance.  Again,  a  purely 
one-dimensional  solution  was  generated.  Oscfffatfons  matched  a  first  longitudinai  mode  with 
overtones.  Next,  the  liquid  was  radially  put  into  the  first  centimeter  of  distance.  This  is  also  a 
one-dimensional  problem,  with  the  profiles  uniform  in  the  axial  direction.  A  first  radial  mode 
was  generated.  In  both  cases,  accurate  solutions  were  obtained  with  a  fairly  coarse  grid. 
Mass  was  conserved  almost  exactly.  There  were  energy  errors  due  to  the  work  term.  The 
implicit  formulation  reduces  these  errors  but  does  not  eliminate  them.  For  a  coarse  grid,  the 
energy  error  was  a  small  fraction  of  1%.  For  the  fine  grids,  the  energy  error  was  near  1%. 
The  reason  the  energy  error  gets  worse  for  a  fine  grid  is  due  to  the  procedure  for  choosing  a 
time  step.  Values  are  generated  every  0.002  ms.  For  a  coarse  grid,  the  time  step  from  the 
Courant  condition  would  be  larger  than  this.  Since  the  code  hits  the  output  times  exactly,  it 
takes  a  smaller  step  than  required  for  stability  and  so  has  improved  accuracy.  For  the  finer 
grids,  with  much  smaller  control  volumes,  the  time  step  would  have  to  be  made  smaller  to  get 
equivalent  accuracy.  A  1%  error  in  the  energy  conservation  is  judged  acceptable. 

To  obtain  an  actual  two-dimensional  problem,  2.5  g  of  liquid  was  axially  and  radially  put 
into  the  first  centimeter.  Figure  6  shows  the  resulting  pressure,  measured  at  the  wall  in  the 
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middle  of  the  chamber.  The  pressure  from  the  corresponding  uniform  distribution  is  shown  for 
comparison.  Figure  7  shows  a  blowup  of  part  of  the  curve  for  three  different  grids.  Grids  with 
higher  resolution  in  the  radial  direction  were  used  since  the  radial  modes  are  more  important 
than  the  longitudinal  modes.  Figure  8  shows  the  Fourier  transform  of  the  pressure  curve 
between  1  and  2  ms  for  the  different  grids.  The  major  peak  is  the  first  radial  (28  kHz),  and 
the  second  highest  peak  is  the  first  longitudinal  (23  kHz).  Since  the  curves  are  not  pure 
sinusoidals,  higher  overtones  are  also  generated.  The  first  radial  becomes  more  dearly 
defined  as  the  grid  is  refined.  However,  even  the  coarse  grid  shows  the  main  features  of  the 
solution. 

To  check  the  injection  algorithm,  the  chamber  is  filled  with  gas  at  100  MPa  and  2,468.5  K. 
Then  5  g  of  liquid  is  injected  over  the  first  millisecond  through  a  circular  opening  1 .0  cm  in 
radius.  Mass  conservation  is  exact  for  all  the  grids  tried.  To  check  energy  conservation,  the 
chemical  energy  of  the  injected  liquid  must  be  added  to  the  internal  energy  of  the  gas  in  the 
chamber.  In  addition,  the  injected  liquid  does  work  on  the  gas  in  the  chamber,  compressing  it 
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into  a  smaller  volume.  When  this  is  taken  into  effect,  the  energy  error  is  much  less  than  1%. 
The  first  radial  and  first  longitudinal  modes  are  set  up  as  before.  However,  the  amplitudes  are 
smaller,  since  there  is  much  less  liquid  accumulation. 

The  primer  model  was  also  checked  by  injecting  5  g  of  hot  gas  over  the  first  millisecond. 

A  very  small  first  longitudinal  mode  is  set  up.  This  is  caused  by  the  injected  gas  displacing 
gas  already  in  the  chamber. 

Some  unexpected  results  also  were  observed.  When  liquid  was  radially  put  into  the  first 
centimeter,  a  clean  first  radial  mode  was  observed.  Next,  a  smaller  amount  of  liquid  was 
radially  put  into  the  first  0.25  cm.  The  liquid  was  again  assumed  to  be  in  lOO-pm-diameter 
droplets.  Figure  9  shows  the  observed  frequencies.  The  first  four  radial  modes  were  excited. 
For  a  fine  grid,  the  main  frequency  is  the  third  radial  mode. 

At  time  zero,  the  liquid  combusts  rapidly  and  generates  a  high  pressure  region.  This 
begins  to  move  the  fluid  away  from  the  centerline.  Eventually,  enough  fluid  has  been 
accelerated  toward  the  wall  that  the  pressure  drops,  and  combustion  slows  drastically. 
However,  because  of  momentum,  the  fluid  continues  to  move,  and  the  pressure  drops  below 
the  initial  pressure.  The  flow  reverses,  and  combustion  picks  up  again.  Because  of  the  high 
burning  rate  of  the  propellant  this  process  can  generate  fairly  large  pressure  oscillations 
before  there  are  any  reflections  from  the  chamber  wall.  The  smaller  the  liquid  region,  the 
higher  the  frequency  of  the  generated  waves.  For  the  present  problem,  the  natural  frequency 
of  the  liquid  region  is  close  to  the  third  radial  mode  of  the  chamber,  and  this  mode  is 
preferentially  pumped.  The  coarse  grid  cannot  resolve  the  steep  pressure  gradient  in  the 
liquid  region.  So  as  the  grid  is  refined,  the  main  effect  is  to  increase  the  resolution  of  the  third 
radial  mode. 

I  also  tried  the  same  initial  setup  but  with  50-(xm-diameter  droplets  (Figure  1 0>.  The 
oscillations  are  larger.  But  the  natural  frequency  of  the  liquid  region  has  been  changed,  and 
now  the  first  radial  mode  is  preferred. 
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Figure  1 0.  Cl 


In  liquid  rockets,  the  propellant  is  usually  injected  as  uniformly  as  is  practical.  In  liquid 
guns,  on  the  other  hand,  the  propellant  is  usually  injected  through  an  annulus  that  is  small 
compared  to  the  chamber  area.  Hence,  higher  frequencies  can  be  generated. 

42  Ooen  Chamber.  An  opening  0.5  cm  in  radius  is  put  into  the  right-hand  wall.  For  a 
first  test,  gas  is  put  into  the  chamber  at  100  MPa  and  2,000  K,  and  the  exit  pressure  just 
outside  the  chamber  is  put  at  90  MPa.  Figure  1 1  shows  the  outflow  rate.  For  comparison,  a 
solution  is  also  given  for  a  lumped  parameter  model  (Coffee  1988).  That  is,  the  conditions  in 
the  chamber  are  assumed  to  be  uniform,  and  the  flow  velocity  out  is  given  by  assuming 
steady-state  isentropic  flow.  For  the  detailed  equations,  there  is  a  time  delay  as  the  flow  gets 
going.  The  flow  does  not  get  as  large  as  for  the  lumped  parameter  model  since  there  is  a 
delay  in  getting  fluid  from  inside  the  chamber  to  the  exit  plane.  This  means  the  chamber 
takes  longer  to  empty.  However,  the  time  scale  is  similar.  The  exit  flow  is  almost  grid 
independent.  There  are  no  noticeable  pressure  oscillations  in  this  case. 

Next,  uniform  injection  at  the  left  wall  was  implemented  at  a  rate  of  5  g/ms.  Figure  12 
shows  the  resulting  outflow  rates.  The  rate  is  almost  grid  independent  The  model  does 
reach  a  steady  state,  where  the  outflow  matches  the  inflow.  There  are  some  very  small 
oscillations.  Pressure  waves  that  hit  the  right  wall  reflect,  but  pressure  waves  that  hit  the 
opening  exit,  and  a  rarefaction  wave  propagates  back  into  the  chamber.  The  two  effects  tend 
to  cancel,  and  it  is  difficult  to  start  up  oscillations. 

The  same  basic  setup  was  created  with  choked  outflow.  The  injection  rate  was  increased 
to  7.25  g/ms  to  compensate  for  the  higher  outflow.  For  this  case,  large  oscillations  were  set 
up.  The  choked  flow  boundary  acts  more  like  a  reflecting  wall.  Figure  13  shows  the  pressure 
at  the  chamber  wall,  and  Figure  14  shows  the  Fourier  transforms  of  the  wall  pressure  for 
different  grids.  There  is  a  longitudinal  mode.  Since  the  outflow  is  only  near  the  centerline, 
conditions  are  not  uniform  in  the  radial  direction.  Once  the  flow  becomes  disturbed,  the  first 
radial  mode  is  preferred.  However,  a  fine  grid  is  necessary  to  pick  up  the  radial  mode.  There 
is  also  a  very  low  frequency  around  1  kHz.  This  is  probably  due  to  bulk  oscillations  of  the 
material  in  the  chamber  as  a  whole. 
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Next,  injection  through  an  opening  1  cm  in  radius  in  the  left  wall  was  considered.  Figure 
15  shows  the  Fourier  transform  for  the  pressure  at  the  middle  wail.  With  the  coarse  10x10 
grid,  there  are  no  oscillations.  However,  the  20  x  20  grid  shows  the  proper  behavior.  The 
pressure  shows  a  clean  first  radial  mode  with  overtones.  That  is,  the  pressure  is  not  a  pure 
sinusoidal  wave,  so  higher  frequencies  appear  at  multiples  of  the  first  radial  frequency.  This 
can  be  distinguished  from  the  higher  radial  modes,  which  have  slightly  different  frequencies. 

Figure  1 6  shows  the  corresponding  result  for  choked  flow.  The  oscillations  are  almost  the 
same  magnitude.  The  coarse  grid  shows  only  the  longitudinal  modes.  The  finer  grids  show 
primarily  the  first  radial  mode  with  overtones. 

Next  the  same  amount  of  liquid  was  injected  through  a  hole  0.25  cm  in  radius.  Figure  17 
shows  the  Fourier  l/ansforms.  For  the  20  x  20  grid,  the  highest  frequency  is  the  first  radial. 
There  is  a  frequency  corresponding  to  twice  the  first  radial,  and  a  slightly  lower  frequency 
corresponding  to  the  second  radial.  For  the  40  x  40  grid,  it  is  primarily  the  second  radial 
mode  that  appears.  The  60  x  60  grid  shows  results  very  close  to  the  40  x  40  grid.  Since  the 
liquid  appears  in  a  small  volume,  it  will  generate  pressure  pulses  that  may  excite  higher 
modes.  It  takes  a  fine  grid  to  resolve  this.  The  20  x  20  case  takes  about  5  min  to  run  on  the 
U.S.  Army  Ballistic  Research  Laboratory  (BRL)  Cray  X-MP,  while  the  60  x  60  grid  takes  over 
7  hr. 

Figure  18  shows  the  early  pressure  at  the  control  volume  on  the  centerline  next  to  the 
injector.  The  time  is  too  early  for  any  reflected  waves  to  return  from  the  walls.  The  large 
oscillations  are  generated  from  the  inertial  confinement  of  the  liquid.  A  fine  grid  can  better 
resolve  the  very  steep  pressure  gradient  in  the  small  volume  of  liquid.  However,  this  behavior 
may  owe  something  to  the  simplifications  in  the  model.  The  model  assumes  that  liquid  is 
injected  at  a  steady  rate,  independent  of  the  pressure  by  the  inlet.  The  liquid  breaks  up 
instantaneously  into  droplets  and  also  ignites  instantaneously.  In  reality,  there  must  be  some 
delay  while  the  liquid  breaks  up  and  ignites.  The  combustion  front  will  start  on  the  outside  of 
the  jet  and  move  in.  The  model  merely  shows  that  inertial  confinement  is  capable  of 
generating  large  pressure  waves.  Whether  this  happens  in  reality  is  still  to  be  determined.  In 
the  futur  more  sophisticated  jet  breakup  models  may  be  implemented  in  the  code. 
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Figure  15.  Chamber  With  Outflow.  Injection  Through  a  1  -cm-Radius  Hole.  Fourier 
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Figure  16.  Chamber  With  Choked  Outflow.  Injection  Through  a  1 -cm-Radius  Hole. 
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12.0-1 


Figure  17.  Chamber  With  Outflow.  Injection  Through  a  0.25-cm-Radius  Hole.  Fourier 
Transform  of  Wall  Pressure.  20  x  20  Grid  (line).  40  x  40  Grid  (dot). 

60  x  60  Grid  (dash). 


igure  18.  Chamber  With  Outflow.  Injection  Through  a  0.25-cm-Radius  Hole.  Pressure  in 
Front  of  Injector.  10x10  Grid  (line).  20  x  20  Grid  (dot!.  40  x  40  Grid  (dash). 
60  x  60  Grid  (dot-dash). 


43 


Figure  1 9  shows  the  corresponding  Fourier  transforms  for  the  choked  outflow.  The  major 
frequency  is  the  first  radial.  It  is  not  clear  why  the  behavior  is  so  different  for  choked  and 
unchoked  flow.  It  is  clear  that  rapid  injection  of  liquid  into  a  small  volume  can  cause 
complicated  behavior. 

5.  TURBULENCE 

With  the  large  velocities  in  the  gun,  flow  is  unlikely  to  be  laminar.  So  a  turbulence  model 
is  included  as  an  option  in  the  code. 

Turbulence  is  a  three-dimensional  phenomenon.  Relatively  large  eddies  are  developed  in 
the  flow.  These  eddies  decay  into  smaller  eddies.  At  a  very  small  length  scale,  the  physical 
viscosity  turns  the  eddies  into  heat  Normally,  the  physical  viscosity  is  not  the  rate-limiting 
step,  so  the  value  of  the  physical  viscosity  is  immaterial. 

The  k-e  model  of  turbulence  is  used  in  the  code  (Jones  and  Launder  1972;  Ng  and 
Spalding  1972;  Hinze  1975;  Bradshaw  1978;  Kollmann  1980;  Jones  and  Whitelaw  1982).  The 
simpler  mixing  length  models  do  not  take  into  account  the  history  of  the  flow,  which  is 
important  for  the  transient  problems  of  interest.  The  more  complicated  Reynold’s  stress 
models  are  not  really  developed  enough  for  practical  use.  All  the  turbulence  models  are  only 
approximations.  The  k-e  model  is  normally  adequate  for  internal  flows.  It  does,  however, 
ignore  some  effects,  such  as  curved  streamlines,  it  is  also  very  poor  at  predicting  flow 
separation. 

The  effect  of  turbulence  is  considered  in  the  momentum  equations.  There  is  also  a 
turbulent  transfer  of  energy.  However,  for  the  problems  of  interest  the  temperature  does  not 
vary  much,  and  this  effect  should  be  unimportant. 

Because  the  full  three-dimensional  equations  are  required,  tensor  notation  is  used.  That 
is,  repeated  subscripts  indicate  summation.  The  subscript  1  indicates  the  x  direction,  which  is 
assumed  to  be  the  main  direction  of  flow  (if  there  is  one).  The  subscript  2  indicates  the  y 
direction,  which  is  assumed  to  be  the  main  direction  of  shear  (if  there  is  one).  The  subscript  3 
indicates  the  z  direction. 
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Figure  19.  Chamber  with  Exit  Hole  -  Choked.  Injection  Through  Small  Hole. 

20  x  20  Grid  (Line!.  40  x  40  Grid  (Pott.  60  x  60  Grid  (Dash). 

The  vast  majority  of  work  has  been  for  incompressible  flow.  First,  the  equations  for  this 
case  are  derived.  Then  the  equations  are  generalized  to  compressible  flow. 

5.1  Incompressible  Flow.  The  flow  Is  assumed  to  be  incompressible  and  the  viscosity 
constant  The  continuity  equation  reduces  to 


dVj/dXj  «  0. 


The  momentum  equations  are  given  by 


9v,/9f  +  dV'Vj/dXj  m  - ( 1  /p  ){dp/dx,  +  dXj,/dXj}, 


with 


iji  »  -n(9v,/3xy  +  9vy/3x()  +  2/3p.8;/9v*/9x* 

■  -Viid^/dXj  +  dVj/dx,) . 


Using  the  kinematic  viscosity  (cm2/s), 


v  -  p/p, 

the  momentum  equations  can  be  rewritten  as 

dv,/dt  +  dV'Vj/dXj  ■  -(  1  /p)dp/dx, 

+  d/dXjfrldvJdXj  +  dVjIdx,]}  .  (7) 

The  basic  idea  is  to  split  the  flow  into  a  mean  part  and  a  fluctuating  part  Since  transient 
problems  are  of  interest,  an  ensemble  average  is  used.  That  is,  the  problem  is  repeated  N 
times,  where  N  is  a  large  number,  and  the  mean  is  defined  by 

7(f)  -  Zf(t)/N. 

For  any  given  run  of  the  problem,  /  -  7  +  V  .  where  V  is  the  fluctuating  component.  By 
definition,  the  mean  of  the  fluctuating  part  is  zero. 

Now  the  mean  of  the  momentum  equations  is  taken.  The  only  term  that  is  nontrivial  is 

v^/j  •  {v,  *  v,' )  [Vj  *  v/ ) 

-W  +  v/vj  +  v,v/ 

■^♦77- 

Note  that  a  mean  quantity  can  be  moved  out  from  under  an  overbar  and  that  the  mean  of  a 
fluctuating  component  is  zero. 
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The  momentum  equations  become 

dv,/dt  +  3 vyj/dXj  =  -( 1  /p )dp/dx,  +  3/3x;{v[3^ /dXj  +  3v^/3x,]  -  v/  Vj' } .  (8) 


The  new  term,  the  product  of  the  mean  of  the  fluctuating  velocities,  is  called  the  Reynold’s 
stress  tensor.  This  is  the  only  term  by  which  the  turbulent  fluctuating  velocities  affect  the 
mean  flow. 

The  continuity  equation  is  simply 

3 Vj/dXj  ■  dVj/dXj  ■  3v//3 Xj  -  0. 

For  the  k-e  model,  the  quantity  of  interest  is  the  kinetic  energy  of  the  turbulent  fluctuations, 
defined  by 


k  m  0.5  v/  v/  • 

Note  that  all  information  about  the  directional  intensity  of  the  turbulence  is  lost.  Basically  the 
turbulence  is  assumed  to  be  isotropic  (the  same  in  all  directions).  This  is  expected  to  be  true 
at  small  distance  scales,  where  the  turbulence  is  dissipated.  However,  the  turbulence  is 
created  at  large  distance  scales  and  often  in  a  preferred  direction.  So  information  that  is 
important  for  some  problems  is  lost. 

To  obtain  a  turbulent  kinetic  energy  equation,  multiply  the  instantaneous  momentum 
Equation  7  by  v/  and  take  the  mean.  The  time-dependent  term  becomes 


«r/av>/af-  v/rf/dt  +  v/dv//Bt 

■  0  +  0.5  0v,//  v/  *  Bk/dt. 


The  advection  term  becomes 
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where 


V,'  dv,Vj/dXj  ■  v/dv,v//dxt  *  v/ dv/  Vj/dXj  +  V,'  dv/  v/ /dXj> 


1st  part  -  v/  v/dv/dx,  *  v/  v,d v/ /dx,  -  v/  v/  dv,/d xJf 


2nd  part  -  v/  v/  dv/dx,  +  v/  Vjdv/ /dx,  "  v/  Vjdv/ /dx, 

•  Vjdk/dXj  -  dVjkldXj  -  kdv/dx,  •  dVjk/dx ), 


3rd  part  -  v/  v/dv/  /dx,  ♦  v/  v/dv/  /dx, 

*  0  +  0.5  v/  d  v,'  v,  /dXj 
■  0.5  dv/  v,  v,  !<*x,  -  0.5  v/  v/  dv/ /dXj 
-  0.5  dv/v/~v//dxt. 

The  pressure  term  becomes 

-  (  v/  /p )  a  p/dx,  »  -(1/p)  {  dv,'  p/  dx,  -  pdv/ /dx,  } 

-  -( Mp)dv/’p/dx ,. 

The  viscosity  term  becomes 

v/  a/axy{  v[av,/axy  +  dv/ /dx,  +  dv,/dx,  +  dv/ /dx,]  }«vr/  d2v//dx/ 
+  V  V,'  d3v/ /dXjdx,  «  V  V,'  d2v/  /dx/  +  0  -  vd'k/dx/  -  v{dv/ /dx/2. 


Putting  it  all  together, 


(9) 


dk/dt  +  dVjk/dXj  -  -  v,'  Vj' d  V,ldXj  -  d/dXj{Vj'  v/  v/  /  2 

+  SyPV,'  /p  -  vdk/dXj)  -  v(9v/  /dXj)* . 


The  left  side  of  the  equation  is  now  standard.  There  is  a  time-dependent  term  and  an 
advection  term.  The  right  side  consists  of  a  production  term  (turbulence  is  created  from  the 
mean  flow),  three  diffusion  terms,  and  a  dissipation  term.  The  last  term  is  given  the  notation 


e  -  v(dv,'  IdXj)2. 


Normally,  the  physical  viscosity  terms  are  negligibly  small.  However,  for  this  last  term, 
gradients  in  the  fluctuating  velocities  are  considered,  and  this  has  a  small  length  scale.  This 
term  is  the  only  term  that  destroys  turbulent  kinetic  energy. 

Older  theories  are  based  on  the  idea  of  a  mixing  length,  which  is  a  length  scale  for  the 
turbulent  eddies.  The  relation  is 

1  -  k'*U  . 


Under  the  k-e  theory,  the  mixing  length  can  change  with  time  and/or  space,  rather  than  being 
held  constant 

The  problem  now  is  to  find  expressions  for  the  various  correlations  in  Equation  9.  So 
consider  a  very  simple  test  case.  Assume  a  steady-state  shear  flow.  That  is,  the  only 
nonzero  mean  velocity  is  in  the  x  direction,  and  the  only  nonzero  space  derivative  of  the  mean 
velocities  is  in  the  y  direction.  Also  assume  that  the  diffusion  term  is  negligible.  Then 
Equation  9  reduces  to 


Define 


-  v,'  Vj  dv,/dx2  ■  e. 


*  (  V,  v2  Ik)2. 
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Then 

e2  ■  cuk2{dvjdx2)2 , 

or 

e  -  cu{k2le)(d7t/dxst)2, 
and 

-  v/  v-/  -  cu(k2/t)  dvjdx2. 

Define  the  turbulent  viscosity 

v,  -  cu{k2/e). 

Then 

-  v,V  •v,dfi/dx2, 

or  after  much  effort,  the  Reynold’s  stress  now  looks  just  like  the  molecular  stress  term,  only 
with  a  turbulent  viscosity  rather  than  a  physical  viscosity. 

So  now  the  major  assumption  is  made  that  in  general  flows,  the  Reynold’s  stress  looks 
like  a  viscosity  term.  That  is, 

-  v,'  v/  -  v,[3  7,/dXj  +  d  Vj/dx,]  -  (2/3)5tjk. 

The  last  term  is  required  to  balance  the  equation.  That  is,  if  i = /,  then 

-  v,'  v,'  »  -2  k. 


The  molecular  stress  also  breaks  into  the  following  two  parts:  the  inhomogeneous  term 
(dependent  on  direction)  which  becomes  the  viscosity  term,  and  the  homogeneous  term 
(independent  of  direction)  which  is  called  the  pressure. 

Now  the  Reynold’s  stress  formula  is  substituted  into  the  mean  momentum  Equation  8,  and 
the  result  is 
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dv,/dt  +  dv,vj/dXj  ■  -(1  / p)3p/3x ,  -  (2J3)dk/dx, 

+  3/3xy{  (v+v,)  [3^/3xy  +  3i^/3x,]  }.  (10) 

So  there  is  a  turbulent  pressure  term  and  a  turbulent  viscosity  term.  The  turbulent  pressure  is 
normally  small  compared  to  the  actual  pressure.  In  most  cases,  the  turbulent  viscosity  is 
much  larger  than  the  physical  viscosity.  One  exception  is  in  a  boundary  layer,  such  as  near  a 
wall.  The  flow  will  have  a  laminar  sublayer,  as  well  as  a  region  where  laminar  and  turbulent 
effects  are  about  the  same  size  (see  below). 

For  the  turbulent  kinetic  energy  Equation  9,  first  consider  the  production  term.  This  now 
becomes 


Pk  ■  -v,'  v/  3  7,/dXj 

-  v13^/3xy[3^’/3xy  +  3^/3x)]  -  (2/3)  kd  v,/dx, 


or 

Pk  *  v,3  v^/3xy[3  v^/3xy  +  3^/3x,]. 

The  diffusion  term  is  very  messy  and  is  expected  to  be  of  secondary  importance.  So  the 
assumption  is  made  that  diffusion  works  like  molecular  diffusion.  The  result  is  the  equation 

dk/dt  ♦  3  Vjk/dXj  »  3/3xy[  (v,/o*)  3k/3xy]  +  Pk  -  e,  (11) 

where  a„  is  a  constant  to  be  determined  from  experimental  data. 

An  equation  for  the  dissipation  is  still  required.  This  can  be  derived  using  the  same 
general  procedure  as  for  the  turbulent  energy  equation.  However,  the  resulting  equation  has 
messy  correlations  that  cannot  easily  be  simplified.  Instead,  the  standard  procedure  is  to  just 
apply  dimensional  analysis  and  get  an  equation  following  the  kinetic  energy  Equation  1 1 . 

That  is, 
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(12) 


9e/3 1  +  3  vje/dXj  «  3/3xy[  ( vt/cf)3e/3xy] 

*  cti(t/k)Pk  -  ct2(e2/k). 


The  constants  are  determined  by  looking  at  expe.imentai  data  for  some  simple  cases.  These 
cases  are  normally  homogeneous,  decaying  turbulence  behind  a  grid,  near  wall  turbulence  in 
local  equilibrium,  and  free  shear  flows.  The  standard  set  of  values  is 

0.09  1.44  1.92  1.0  1.3. 


The  above  standard  model  generally  does  well  for  internal  flows,  as  long  as  body  forces  are 
negligible  and  the  streamlines  are  not  curved.  However,  some  people  believe  that  it  is 
fundamentally  wrong  and  only  works  because  the  constants  have  been  adjusted  to  match 
experimental  data.  Unfortunately,  this  makes  it  a  little  risky  to  extrapolate  to  gun  conditions, 
since  almost  all  measurements  have  been  made  at  low  pressures. 

5.2  Compressible  Flow.  The  case  of  actual  interest  is  compressible  flow.  The  kinematic 
viscosity  is  assumed  to  be  the  same  as  the  mean  kinematic  viscosity.  The  continuity  equation 
is  now 


3p/3f  +  dpVj/dXj  -  0.  (12) 

The  momentum  equations  are  given  by 

3pv,/3 1  *  dp^Vj/dXj  ■  -dp/dx,  -  dTj'/dXj,  (13) 

with 

Xj,  ■  -p(3v,/3xy  +  dVj/dx,)  +  2J3\iB,jdvk/dxk 
■  -pv(3v,/3xy  +  3vy/3x,)  +  2/3pvS(y3v)t/3xilI. 


52 


Note  that  the  last  term  of  the  stress  tensor  is  no  longer  zero.  The  standard  ensemble  average 
introduced  above  is  used  for  the  density. 


p  -  p  +  p;. 

However,  for  the  velocity,  a  mass-averaged  quantity  is  introduced— that  is, 


where 


or 


This  implies  that 


Then 


v,  •  v,  +  v/  -  v,  +  v,"  , 


V,  -  p  v,/p 


p*/  -  Pvr 


P  V'  «  p  V,  -  p  v,  -  0. 


p  w,  -  p^ +  p  yi"  +  p  v>"  22 ♦  pi'/'  vj" 

* 

-  P v, *  v,pVi"  +  ^p  V,"  +  py/;  v/' 

■py.ypy/7  V'. 

Take  the  mean  of  the  continuity  equation, 


dp/9f  +  9p  Vj/dXj  -  0 . 
Take  the  mean  of  the  momentum  equations, 


dp  v,/3f  +  dp  ViVj/dxj  -  -dp/dx,  +  d/dx^pv^d  v,/dXj  +  d  Vj/dx,  j 

-  2/3 p  v  8f/9  v„/Bxh  -  p  v,"  v )"  j . 
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Now  define  the  turbulent  kinetic  energy  as 

k  -  0.5  v/V'. 

As  before,  assume  that  the  Reynold's  stress  can  be  written  as 

-pv/'v/'  ■  pv,l9  v,/9xy  +  9  vy/9x(] 

-  (2/3)5<yp  v,9  -  (2/3)8,yp/c. 

Then  the  mean  momentum  equations  become 

dpv,/dt  +  dp  v,vy/9xy  ■  -9p/9x,  -  (2/3)9pfc/9, 

+  9/9xy{p(v+vf)[9  v,/9xy  +  9vy/9x(] 

-  ^ptv+vJS^v^/ax,}.  (14) 


Next  take  the  instantaneous  momentum  Equations  13,  rearrange  using  the  continuity 
Equation  12,  and  plug  in  the  stress  term  to  obtain 

p9v,/9f  +’p  VjdV'/dXj  -  -9p/9x,  +  9/9xy(pv[9vi/9xy  +  9  vy/3x(] 

-  2/3pv,y9v*/9x*). 

Multiply  by  v,"  and  take  the  average.  The  time-dependent  term  becomes 

p  V,"  9  V,tdt  »  p  V,"  9v,"/9f  »  pdk/dt. 

The  advection  term  becomes 


54 


p  V,"  VjdV/ldXj  »  p  V,"  Vj"  dv,/dXj  *  p  V,"  Vjdv," /dXj 
+  pv,"  Vj"  dv/' /dXj 

■  P  V,"  v"d  V,/dx{  +  p  Vjdk/dxJ  +  P  0.5  v"  dv"  V," /dxr 

The  pressure  term  becomes 

-  V,"  dp/dxr 

The  first  part  of  the  viscosity  term  becomes 

v,"  d/dXjipv  JdvJdXj  +  dVj/dx,) }  ■  3/9x;{pv  v,"  {dv," /Bx,  +  dv"  Id x,) ) 

-  p v{dv," /dXj  +  dv/' /dx,)dv," /dxt. 

The  second  part  of  the  viscosity  term  becomes 

-  v,"  d/dXj{2/3pvdlldvk/dxk}  *  -  v,“  d/dx,{2!2  p vd vjdxk) 

-  -dldx,{2IZpvv,"  dv"  !dxk)  -  2/3pvdv" /dxkdv," /dxr 


Combining 
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p9fc/9f  +  P  Vjdk/dXj  -  -p  W  v"  3  V,/d Xj 

~  p0.5  v/'dv,"  v,"  IdX]  -  V,"  9p/9x, 

+  9/9Xy{pV  V,"  (3  V,"  IdXj  +  3Vy"/3x,)} 

-  d/dX'jMpvv/'dv" /dxj 

-  pv(3v,"/3Xy  *  dv/'/dx, )**,'"* 

-  2/3f>vdvk" /dxkdv," /dx,.  (15) 


This  is  somewhat  messier  than  Equation  9  for  incompressible  flow,  but  the  basic  structure  is 
the  same.  The  first  line  gives  the  time-dependent  term,  the  advection  term,  and  the 
production  term.  The  next  three  lines  are  diffusion  terms.  Again,  for  lack  of  a  better  idea,  it  is 
assumed  that  diffusion  works  like  molecular  diffusion.  The  last  two  lines  are  dissipation  terms, 
and  ail  of  them  are  set  equal  to  pe.  So  the  result  is 


and 


P„  »  v,[9  v^dxj  +  9  VjIdX'ldv./dXj 

-  (2/3)vtdvk/dxkdvk/dxk  -  (2/3 )  kdvk/dxk 

e  -  -v(3v// /9xy  +  dv," /dx,)dVl" /dx, 

-  2/3v9v*"/3x,9v*"/3x*. 


Rearranging  the  left  side  of  Equation  15  using  the  continuity  equations  and  plugging  in  the 
indicated  expressions  on  the  right-hand  side  results  in 

9p/c/9f  +  9p  Vyk/3Xy  »  9/9Xy[(p.f/o„)3k/9Xy]  +  p  Pk  -  pe.  (16) 
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Similarly,  the  dissipation  equation  becomes 


dp e/df  *  dp  Vje/Bxj  «  d/dXy[  ( pf/ot) de/dXy]  +  ct1  (e/ k)p  Pk  -  ct2p(  e2/ k) .  (17) 

Since  there  is  almost  no  data  for  compressible  flows,  the  constants  are  given  the  same  values 
as  for  incompressible  flow. 

5.3  Implementation.  This  new  set  of  equations  must  be  solved.  From  now  on  all 
quantities  are  assumed  to  be  the  appropriate  means,  so  overbars  and  underlines  will  be 
omitted.  Taking  Equation  14  and  writing  it  for  the  x  momentum  gives 

d(p v,)/df-  -  d(pv/)/dx  -  d(pv,vy)/dy  -  dp/dx  -  (2/3)dp/c/dx 
+  d/dx[4/3  ( p  ♦}!,)  d  V,/dx  -  2/3  (p  +  p()dvy/dy] 

*  d/dy[ (p  +  PjJdv^/dy  ♦  (ji  +  pf)dvy/dxj.  (18) 


Applying  the  divergence  theorem, 

d/dfJy(pvJcW«  -fs{pvMvJ  +  pvxvyJ  +  [p  +  ( 2/3 ) p k] /) nPS 

♦  J5([4/3(p  +  p,)dv,/dx  -  2/3(p  +p,)dv//dy]/ 
+  [(fn-|if)3v,/dy  +  (\L  +  \it)dvy/dx]j)ndS. 

Similarly,  the  momentum  equation  in  the  radial  direction  is 

d(pvy)/df  -  -  d(pv,vr)/dx  -  d(pvy vr)/dy  -  dp/dy  -  (2/3)dp/c/dy 
♦  d/dx[(|l  +  |J.f)dvjt/dy  +  (n+pf)d  Vy/dx] 

+  d/dy[4/3(p +p()dvy/dy  -  2/3  (  p  +  p,)  d  vj  dx] . 


(19) 


Applying  the  divergence  theorem. 


3/3fJy(pvy)dV-  -js[pvxvyl+pvyvyJ*[p*(2J2)pk]J)ndS 
*  +  (\l*\L,)dvy/dx]l 

+  [4/3(p+p,)3vy/3y  -  2/3(p  +  p,)3v't/3x]/)/JdS. 


Next  consider  the  turbulent  production  term.  For  a  two-dimensional  problem,  this  becomes 

pP,  -  p,[(4/3)(3v,/3x)2  ♦  (3 v-,/3y)2  +  (3vy/3x)2 
♦  (4/3)(3vy/3y)2  +  2(3vjr/3y)(3vy/3x) 

-(4/3)  (dvx/dx)(dvy/<}y)  ] 

-  (2/3)p*(3v,/3x  +  3vy/3y). 

The  turbulent  energy  Equation  16  becomes 


or 


dpk/dt  »  -3p  vxk/dx  -  dpvyk/dy  +  3/3x(  pt3/c/3x) 

+  3/3y(pf3/c/3y)  ♦  p  Pk  -  pe 

d/dtfv(pk)dV~  fv{pPk-pt)dV-  fs(pvxkl  +  pvykj)ndS 
+  J*g(p,3/c/3x/  +  p73fr/3y/)  ndS. 


The  dissipation  Equation  17  becomes 

3pe/3f  »  -3pv,e/3x  -  3pvye/3y  +  3/3x[(pf/1.3)3e/3x] 

♦  3/3y[(p,/1.3)3e/3y]  +  1.44(e/fr)p  Pk  -  1-92 p  (  e2/k) 
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or 


a/arJ^(pe)dV'-  J^[1.44(e//c)pP*  -  1 .92  p  ( e2  /  Ac)  }dV 
-  Js(pv,e/  +  p  vyeJ)ndS 
+  Js[(m/1.3)ae/3x/  +  (n,/1.3)ae/ay/]ndS. 


The  new  quantities  are  put  on  the  scaiar  grid.  There  is  still  the  problem  of  initial  conditions.  A 
gun  firing  cycle  starts  with  stagnant  gas,  so  there  should  be  no  turbulence.  However,  the 
turbulent  terms  cannot  be  zero  because  then  the  turbulent  viscosity  is  undefined.  If  the 
turbulent  kinetic  energy  is  set  to  zero,  but  the  dissipation  is  nonzero,  the  kinetic  energy  term 
will  become  negative.  So  for  initial  conditions,  pk  is  set  equal  to  1 0'2  and  pe  is  set  equal  to 
10"5.  The  density  times  the  turbulent  quantities  is  used  since  the  equations  are  being  solved 
in  conservation  form.  This  makes  the  turbulent  kinetic  energy,  dissipation,  and  turbulent 
viscosity  all  small. 

5.4  Boundary  Conditions.  For  flow  close  enough  to  a  wall,  the  flow  must  be  laminar,  and 
the  turbulent  kinetic  energy  will  be  zero.  However,  while  kinetic  energy  cannot  be  generated 
at  a  wall,  it  can  diffuse  in  and  be  dissipated  at  the  wall.  Hence,  a  common  boundary  condition 
for  the  dissipation  is  to  assume  a  zero  gradient.  That  is,  the  dissipation  in  a  scalar  control 
volume  next  to  the  boundary  equals  the  dissipation  in  the  next  scaiar  control  volume  away 
from  the  wall.  The  turbulent  kinetic  energy  is  then  still  defined  (equals  zero).  At  the 
centerline,  reflection  boundary  conditions  are  used,  and  at  the  exit  from  the  chamber, 
extrapolation  conditions  are  used. 

In  practice,  the  boundary  layers  are  too  small  to  be  resolved  with  a  reasonable  grid,  and 
more  complicated  procedures  are  required.  So  consider  flow  near  a  wall.  Assume  that  the 
radial  flow  (perpendicular  to  the  wall)  is  negligibly  small.  Assume  that  the  axial  gradients  are 
negligibly  small.  These  assumptions  will  be  good  if  we  go  close  enough  to  the  wall,  where 
close  enough  depends  on  the  conditions  of  the  flow.  Under  these  assumptions,  the  only 
nonzero  stress  is 

-  ~(\i  +  \Lt)dvx/dy. 
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Very  near  the  wall,  the  flow  becomes  laminar.  So  the  stress  at  the  wall  is  given  by 


-  -\idvx/dy. 

Assume  that  the  stress  is  constant  through  the  boundary  layer.  Define  the  wall  friction 
velocity  as 

W*  -  J-Xj p  . 


In  this  derivation,  to  be  consistent  with  the  usual  notation,  y  is  the  distance  from  the  wall 
rather  than  the  distance  from  the  centerline.  In  the  laminar  sublayer, 

dvM/dy  -  -xJil 

-  ~(xw/\i)y 
vx  -  pu,2y/n 


Defining 


Vg  «  vx! um  y *  «  t/*y/v 


then  in  the  laminar  sublayer, 

v/  -  y  *  .  (20) 

Next  consider  flow  in  the  turbulent  layer,  where  the  laminar  viscosity  is  negligible 
compared  to  the  turbulent  viscosity.  Then 

3v,/9y-  -  VlV 
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Assume  that  the  kinematic  turbulent  viscosity  is  given  by 

vf  -  Kify, 

where  k  is  a  constant  to  be  determined.  The  basic  idea  is  that  the  turbulent  eddy  size  is 
constrained  by  the  distance  to  the  wall.  The  formula  is  in  good  agreement  with  the 
experiment.  Then 

ku* ydvjdy  -  -xjp  *  u*2 
3v,/9y  »  u’/<y 

or 

<y*dv*/dy*  *  1 
v*  ■  { 1  Ik)  Iny  *  +  C. 

From  experimental  data,  the  usual  constants  chosen  are 

v;  -  2.5 Iny*  *  5  .  (21) 

There  will  be  a  regime  in  between  the  above  regions  where  both  laminar  and  turbulent 
viscosity  is  important.  Various  models  have  been  proposed  for  this  region.  However,  for 
simplicity,  this  will  be  ignored.  Instead,  the  boundary  layer  profile  will  be  assumed  to  change 
instantaneously  from  the  laminar  sublayer  to  the  turbulent  layer.  This  will  occur  when  the  two 
velocities  match,  that  is, 

y*  *  2.5  Iny *  +  5 

or 

y*  =  10.993  . 
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To  implement  this  formulation,  consider  the  first  vector  grid  point  next  to  the  top  wall,  at  a 
distance  y„  from  the  wall.  This  point  is  assumed  to  be  in  the  turbulent  layer,  and  Equation  21 
is  solved  iteratively  for  u. 

Given  the  wall  friction  velocity,  the  distance  from  the  wall  to  the  point  where  the  transition 
from  laminar  to  turbulent  flow  takes  place,  denoted  by  ye,  can  be  computed.  Suppose  that  ye 
is  greater  than  yb.  This  implies  that  the  original  assumption  was  wrong  and  that,  in  fact,  the 
vector  grid  point  is  in  the  laminar  sublayer.  In  this  case,  the  previous  boundary  conditions  are 
used. 

Next  suppose  that  ye  is  much  less  than  yb.  So  the  laminar  sublayer  is  negligibly  small. 
Either  Equation  20  or  Equation  21  can  be  solved  to  obtain  the  velocity  at  yc,  denoted  as  vc. 
This  can  be  taken  as  the  velocity  at  the  wall,  neglecting  the  very  thin  laminar  sublayer. 

Since  a  simulation  normally  starts  with  stagnant  gas,  there  will  be  a  transition  from  the 
laminar  to  the  turbulent  regime.  To  make  the  transition  smooth,  the  wall  velocity  is  actually 
defined  as 


^ -  ve(y,,-ye)/yb  o <ye<yb. 

Boundary  conditions  on  the  turbulent  quantities  are  also  required.  Various  assumptions 
are  possible.  The  simplest  of  the  standard  assumptions  is  that  local  equilibrium  exists  in  the 
turbulent  boundary  layer,  that  is 


e  -  Pk  «  v,(dvx/dy)2 . 

This  means  that  the  production  of  turbulent  kinetic  energy  exactly  matches  the  dissipation. 
The  production  term  simplifies  to  the  right-hand  expression  in  a  boundary  layer.  Using  the 
expressions  derived  above  for  the  turbulent  part  of  the  boundary  layer,  the  above  equation 
becomes 


£  *  KWy(u’/Ky)2  *2.5 u’3/ y. 
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By  the  definition  of  the  turbulent  viscosity, 


or 


v,  ■  0.09  ( Ac 2/e ) 

/r«  u*2/ 0.3. 


Now  consider  the  laminar  sublayer.  By  definition,  k  =  0.  The  simplest  standard 
assumption  is  that  the  dissipation  is  constant  in  the  laminar  sublayer,  or 

£  -  2.5  u*3/ye . 

Now  what  is  required  is  the  integral  over  the  volume  of  the  density  times  the  turbulence 
quantities.  Consider  first  the  top  boundary.  Let  R,  be  the  radius  of  the  cylinder.  Let  Rb  be  the 
radius  up  to  the  bottom  of  the  scalar  control  volume  at  the  boundary,  and  let  Rc  be  the  radius 
up  to  the  crossover  point  between  the  turbulent  and  laminar  regions.  Then  the  control  volume 
can  be  split  into  a  turbulent  and  laminar  section 


Vt  -  dxx(R?  -  R2) 

V;  -  dxn(R,2  -  R2)  . 

The  turbulent  kinetic  energy  is  zero  in  the  laminar  part  and  constant  in  the  turbulent  part.  So 

jypkdV  -  (p  u*2/0.3)  Vr 


The  dissipation  is  constant  over  the  laminar  sublayer,  so 


fvPedV,  -  (2.5p  um3/ye)  V, . 
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Over  the  turbulent  layer 


fvptdV,  -  jv(2.5pu'3/y)dV, 

-  fy‘(2.Spu«/y)2n(Rt-y)dxdy 
*  y» 

-  5 up  u'*dx[R,lny-y]yy\ 

-  5rtp u'3dx[Rt(lny„-lnye)  -  (y„~ye)]. 

The  two  terms  are  added. 

It  is  slightly  simpler  to  work  directly  with  the  turbulent  viscosity.  The  turbulent  viscosity  in 
the  laminar  sublayer  is  zero.  The  turbulent  viscosity  in  the  turbulent  layer  is  given  by  the 
formula 

\it  -  Kp  u * y . 

So 

jv0Apu9ydV 

■  OAp  u*  y2ic(R,-y)dx  dy 

-  Q.Qndxu*  [R,y*/2-y3/3  ]£ 

■  Q.Bitpdxu*[R,{y? -y?)/2  -  (y63-yc3)/3]. 


This  is  divided  by  the  control  volume  to  obtain  the  average  turbulent  viscosity,  and  the 
average  dissipation  is  then  found  from  the  standard  relation 


vt  -  0.09 ( k2/e) . 


The  dissipation  integrals  then  do  not  have  to  be  evaluated. 

Similar  but  slightly  simpler  expressions  can  be  derived  for  the  left-hand  and  right-hand 
boundaries. 
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5.5  Solution  Procedure.  The  two  new  turbulence  equations  are  solved  along  with  the 
previous  equations.  The  first  step  is  still  the  Lagrangian  step.  The  pressure  and  viscosity 
terms  for  the  momentum  equations  are  evaluated,  and  the  velocities  are  updated.  Then  the 
turbulent  energy  production  term  is  evaluated.  The  diffusion  of  the  kinetic  energy  and  the 
dissipation  can  then  be  computed,  and  the  turbulence  terms  are  all  updated.  The  boundary 
conditions  are  also  updated. 

Then  the  Lagrangian  step  is  taken.  The  combustion  term  is  computed.  It  is  assumed  that 
the  turbulence  quantities  do  not  change  during  this  stage. 

Finally,  the  convection  terms  are  computed.  Convection  for  the  turbulence  quantities  is 
done  exactly  the  same  as  for  the  other  scalar  quantities.  The  boundary  conditions  are 
applied,  and  ail  the  quantities  are  updated. 

5.6  Test  Problems.  The  laminar  flow  test  problems  were  now  done  assuming  turbulent 
flow.  Since  there  is  now  a  boundary  layer,  none  of  the  test  problems  reduces  to  a 
one-dimensional  problem.  So  the  two-dimensional  problems  are  the  ones  of  most  interest. 

For  the  closed  chamber  problems,  the  additional  turbulence  had  only  a  minor  effect. 
Qualitatively,  the  results  were  the  same. 

For  the  open  chamber,  the  results  were  more  interesting.  Consider  the  chamber  with 
uniform  injection.  Figure  20  shows  the  new  exit  flow  (comparable  to  Figure  12).  The  outflow 
rate  is  no  longer  grid  independent  because  of  the  boundary  layer  at  the  exit  comer.  However, 
convergence  is  good  with  the  20  x  20  grid.  Small  oscillations  do  eventually  start  in  this  case. 

Next,  the  problem  was  run  with  choked  outflow.  Figure  21  shows  the  Fourier  transforms 
from  4  to  5  ms  (comparable  to  Figure  14).  The  first  radial  mode  now  shows  up  even  with  the 
coarse  grid.  This  is  probably  due  to  the  boundary  layers  disturbing  the  flow.  Because  of  the 
turbulent  viscosity,  the  oscillations  are  not  quite  as  large.  However,  this  also  means  that  the 
coarser  grids  are  a  better  approximation.  The  turbulent  viscosity  will  smooth  out  some  of  the 
very  fine  structure. 
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The  chamber  was  run  with  injection  through  a  large  hole  (Figures  22  and  15).  The  results 
are  similar  to  the  laminar  case  except  that  the  first  radial  does  appear  even  with  the  coarse 
grid  and  the  magnitudes  are  slightly  smaller.  The  same  behavior  occurs  if  the  outflow  is 
choked  (Figures  23  and  16). 

Next  injection  through  the  small  hole  was  considered  (Figures  24  and  17).  The  second 
radial  mode  is  the  primary  frequency.  For  the  40  x  40  grid,  the  oscillations  are  very  small. 

For  choked  flow  (Figures  25  and  19),  convergence  is  very  bad.  The  agreement  is  actually 
better  at  earlier  times.  For  the  finer  grids,  the  oscillations  start  to  die  out  at  the  later  times. 

As  for  laminar  flow,  rapid  injection  into  a  small  volume  can  cause  complicated  behavior. 

6.  GUN  TUBE 

The  model  for  the  gun  tube  is  set  up  almost  exactly  like  the  chamber  model.  The  gun 
tube  is  a  cylinder  with  a  projectile  for  the  right-hand  wall.  The  projectile  acceleration  depends 
on  the  average  pressure  on  its  base.  The  velocity  and  displacement  are  found  by  using  the 
usual  equations  of  motion.  The  grid  in  the  tube  is  attached  to  the  base  and  stretches 
uniformly  as  the  projectile  moves. 

Because  the  gun  tube  volume  will  increase  dramatically  over  the  firing  cycle,  additional 
grid  points  can  be  added  as  the  projectile  moves.  When  the  grid  spacing  becomes  larger 
than  a  specified  value,  a  grid  point  is  added  at  the  end  of  the  time  step.  A  new  scalar  control  ? 

volume  will  normally  contain  parts  of  two  old  scalar  control  volumes.  The  appropriate  fractions 
of  the  old  scalar  quantities  are  combined  to  obtain  the  new  scalar  quantity,  so  the  regridding  is 
conservative.  The  new  velocities  are  found  by  interpolation. 

The  only  new  boundary  condition  is  at  the  interface  with  the  chamber.  The  values  on  the 
gun  tube  grid  just  to  the  left  of  the  tube  are  giver.  e  corresponding  chamber  values.  Some 
adjustment  may  be  necessary,  since  the  grids  in  the  chamber  and  tube  normally  have 
different  spacing.  Similarly,  the  values  on  the  chamber  grid  just  to  the  right  of  the  chamber 
are  given  the  corresponding  tube  values. 
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Frequency  (kHz) 

Figure  22.  Chamber  With  Outflow.  Injection  Through  a  1-cm-Radius  Hole.  Turbulence. 

Fourier  Transform  of  Wall  Pressure.  10x10  Grid  (line).  20  x  20  Grid  (dot 
40  x  40  Grid  (dash 


0.0  20.0  40.0  60.0  80.0  100.0 

Frequency  (kHz) 

Figure  23.  Chamber  With  Choked  Outflow.  Injection  Through  a  1-cm-Radius  Hole. 
Turbulence.  Fourier  Transform  of  Wall  Pressure.  10x10  Grid  (line 
20  x  20  Grid  (dot).  40  x  40  Grid  (dash 
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Figure  25.  Chamber  With  Choked  Outflow.  Injection  Through  a  0.25-cm-Radius  Hole. 

Turbulence.  Fourier  Transform  of  Wall  Pressure.  20  x  20  Grid  (line).  40  x  40 


To  take  a  time  step,  the  pressure  and  viscosity  terms  are  computed  in  the  chamber  and 
then  in  the  tube.  The  boundary  conditions  are  updated.  Then  the  combustion  terms  are 
computed  in  the  chamber  and  tube.  Note  that  droplets  can  enter  into  the  tube  and  combust, 
although  the  amount  of  liquid  in  the  tube  is  generally  small.  Again  the  boundary  conditions 
are  updated.  Finally,  the  convection  terms  in  the  chamber  and  then  the  tube  are  computed, 
the  new  pressures  are  computed,  and  the  boundary  conditions  are  updated.  If  necessary,  a 
new  grid  point  is  added  to  the  tube. 

To  test  the  procedure,  some  simple  test  problems  were  set  up  and  comparisons  made 
with  the  lumped  parameter  gun  code.  The  chamber  was  filled  with  hot  gas,  which  expands 
and  pushes  the  projectile  down  tube.  For  the  first  case,  the  chamber  and  gun  tube  were 
given  the  same  diameter.  For  a  laminar  simulation  (slip  boundary  conditions)  the  problem  did 
reduce  to  one  dimension  (no  radial  variation).  The  chamber  pressure  and  projectile  velocity 
were  identical  for  the  lumped  parameter  code  and  for  the  two-dimensional  simulation. 

Next,  the  chamber  was  made  larger  than  the  gun  tube.  The  chamber  radius  was  3  cm, 
and  the  tube  radius  was  1 .5  cm  (30-mm  gun).  The  gas  was  set  to  an  initial  pressure  of 
100  MPa.  The  integration  was  carried  out  to  6  ms,  at  which  point  the  projectile  was  just  over 
100  cm  down  tube.  The  chamber  grid  was  20  x  20.  At  the  end  of  the  integration,  the  gun 
tube  grid  was  400  x  1 0.  Again,  the  chamber  pressure  and  projectile  velocity  were  identical  to 
the  lumped  parameter  code.  Figure  26  shows  the  gun  tube  velocity  profiles  at  the  end  of  the 
integration.  The  solid  line  is  the  lumped  parameter  code,  using  a  modified  Lagrange  pressure 
distribution.  The  dotted  line  is  flow  along  the  centerline,  and  the  dashed  line  is  flow  along  the 
wall.  The  flow  was  laminar  with  a  slip  boundary  condition.  The  flow  is  more  rapid  along  the 
centerline.  In  fact,  there  is  a  small  recirculation  region  at  the  wail  where  the  flow  enters  the 
tube.  However,  the  flow  does  become  more  one-dimensional  farther  down  tube.  Figure  27 
shows  the  same  curves  for  turbulent  flow.  The  behavior  is  very  different.  The  boundary  layer 
makes  the  flow  much  slower  along  the  tube  wall.  However,  there  is  a  singularity  where  the 
tube  wall  meets  the  projectile,  since  the  gas  at  the  projectile  must  move  axially  at  the  speed  of 
the  projectile.  At  this  corner,  the  fluid  velocity  must  match  the  projectile  speed.  That  means 
that  gas  flows  more  rapidly  down  the  centerline  and  then  moves  radially  up  near  the  projectile 
to  fill  the  space  left  as  the  projectile  moves.  However,  the  pressure  equilibrates  even  though 
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Figure  26.  Gun  Tube  Test  Problem  -  Laminar.  Velocity.  Lumped  Parameter  Model  (line). 
Two-Oimensionai  Model.  Centerline  (dot).  Wall  (dash). 


Figure  27.  Gun  Tube  Test  Problem  -  Turbulent.  Velocity.  Lumped  Parameter  Model  (line). 
Two-Dimensional  Model.  Centerline  (dot).  Wall  (dash). 
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the  velocity  does  not  (Figure  28).  There  is  some  oscillation  near  the  projectile  due  to  the  flow 
disturbance. 


7.  SANDIA  TEST  FIXTURE 

A  simplified  test  fixture  has  been  developed  at  Sandia  National  Laboratories  to  study 
pressure  oscillations.  Figure  29  shows  a  picture  of  the  combustion  chamber  of  the  fixture. 
Liquid  propellant  is  injected  through  the  circular  opening  at  the  left  by  a  differential  piston.  The 
liquid  combusts  and  bums.  Hot  gas  and  possibly  some  liquid  flow  through  the  converging 
orifice  at  the  right  into  atmospheric  pressure.  There  is  a  burst  disc  at  the  end  of  the  injection 
orifice  to  prevent  liquid  propellant  from  leaking  into  the  chamber  prematurely.  There  is  also  a 
burst  disc  at  the  end  of  the  converging  orifice  so  outflow  does  not  start  until  the  combustion  is 
started.  A  primer  initially  prepressurizes  the  combustion  chamber. 

In  the  model,  the  chamber  is  represented  by  a  cylinder.  The  radius  is  taken  to  be  the 
radius  of  the  body  of  the  chamber  (one  inch).  The  length  is  chosen  to  obtain  the  proper 
chamber  volume  (volume  between  the  two  burst  discs).  A  choked  flow  boundary  condition  is 
applied  at  the  right  boundary. 

The  chamber  is  initially  given  a  uniform  pressure  of  22.0  MPa  due  to  the  igniter.  The 
differential  piston  is  driven  by  pressurized  helium.  The  piston  motion  can  be  calculated  from 
the  recorded  helium  pressure.  The  point  at  which  the  burst  disc  breaks  can  be  found  from  the 
sudden  dip  in  the  liquid  pressure.  This  is  taken  as  the  zero  time  in  the  calculation.  In  the 
experiment,  65  cm3  of  propellant  is  injected,  followed  by  water.  From  the  liquid  pressure,  the 
density  of  the  water  and  the  propellant  can  be  calculated.  The  amount  of  liquid  in  the 
reservoir  is  then  calculated,  which  gives  the  amount  of  liquid  injected.  The  derivative  of  the 
amount  of  liquid  injected  gives  the  rate  of  injection  required  by  the  code.  The  code  does  not 
presently  allow  two  different  kinds  of  liquid.  So  at  the  point  where  the  propellant  is  computed 
to  be  exhausted,  the  size  of  the  injected  droplets  is  increased  by  several  orders  of  magnitude. 
This  makes  the  rate  of  combustion  negligible. 
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Figure  29.  Santiia  Test  Fixture  -  Combustion  Chamber. 
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The  right-hand  boundary  is  treated  as  a  solid  wall  before  the  burst  disc  breaks.  When  the 
disc  breaks  (41.4  MPa),  the  boundary  condition  for  the  exit  orifice  becomes  that  of  choked 
outflow. 

Data  were  obtained  for  a  typical  test  (CT-48).  There  are  four  pressure  gauges  in  the 
same  plane  3.5  inches  from  the  injector.  Figure  30  shows  the  recorded  pressures.  All  the 
gauges  show  large  pressure  oscillations,  but  some  gauges  consistently  show  a  higher 
oscillation  level.  Experimentation  has  shown  that  if  the  holes  in  which  the  pressure  gauges 
are  mounted  are  enlarged,  the  magnitude  of  the  oscillations  increases.  So  it  is  felt  that  the 
higher  level  of  oscillations  is  closer  to  what  occurs  in  the  chamber  and  that  the  lower  levels 
are  caused  by  damping  in  the  gauge  hole  (Griffiths  1991). 

Figure  31  shows  the  filtered  chamber  pressures.  There  is  some  variation.  All  of  the 
pressure  gauges  eventually  record  negative  pressures.  Pressure  gauges  do  tend  to  record 
low  when  they  become  hot  (thermal  drift),  so  the  actual  chamber  pressure  may  be  slightly 
higher  than  the  recorded  pressures. 

In  Figure  32,  the  Fourier  transforms  are  given  for  a  1-ms  time  window.  There  is  a  wide 
range  in  the  magnitudes.  However,  ail  of  the  gauges  show  a  fairly  dean  first  radial  mode. 

The  higher  frequendes  are  overtones  due  to  the  fact  that  the  pressure  oscillations  are  not 
pure  sinusoidal  waves.  That  is,  the  higher  frequencies  are  multiples  of  the  first  radial  mode, 
rather  than  the  slightly  different  frequencies  corresponding  to  the  second  and  third  radial 
modes. 

A  model  was  set  up  for  this  fixture.  The  liquid  is  assumed  to  enter  as  200-jim  drops. 
Figure  33  shows  the  filtered  chamber  pressure  from  the  model  measured  at  the  wall  for  three 
different  grids.  The  model  pressure  decreases  substantially  as  the  grid  is  refined.  For  a  finer 
grid,  the  jet  is  less  spread  out,  and  more  of  the  liquid  exits  the  chamber  unburnt.  The 
sensitivity  to  the  grid  size  is  probably  due  in  part  to  the  simplification  in  the  outflow  condition. 
The  actual  fixture  has  a  converging  nozzle,  which  the  model  cannot  presently  duplicate.  If  the 
liquid  is  assumed  to  enter  as  smaller  droplets,  all  the  liquid  burns  in  the  chamber,  and  the 
filtered  pressure  is  not  very  sensitive  to  the  grid  refinement. 
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Figure  31.  Sandia  Test  Fixture.  Filtered  Chamber  Pressures.  P32.  (line).  P33  (dot). 
P34  (dash).  P35  (dot-dash). 
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Figure  32.  Sandia  Test  Fixture.  Fourier  Transforms  -  5  to  6  ms.  P32  (line).  P33  (don. 
P34(dash).  P35  (dot-dash). 


Figure  33.  Sandia  Test  Fixture.  Filtered  Chamber  Pressures.  P35  (line).  Model. 
20  x  20  Grid  (dot).  40  x  40  Grid  (dash).  60  x  60  Grid  (dot-dash). 
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In  Figure  34,  the  unfiltered  pressures  are  given  for  the  three  grids.  For  the  coarse  grid, 
the  agreement  with  experiment  is  reasonable,  although  the  oscillations  are  a  little  too  regular. 
Unfortunately,  for  the  finer  grids,  the  oscillations  decrease  drastically.  This  effect  has  been 
seen  in  some  of  the  test  problems.  When  the  liquid  combusts  in  a  very  small  volume,  it  can 
generate  high  frequency  pressure  oscillations  which  do  not  pump  the  first  radial  mode. 

Figure  35  shows  the  Fourier  transforms  for  the  three  grids.  With  the  finer  grids,  the  higher 
radial  modes  are  excited.  A  number  of  variations  in  the  above  model  were  implemented.  In 
every  case,  when  the  grid  became  fine  enough,  the  oscillations  essentially  disappeared. 

This  problem  may  be  due  in  part  to  the  simplified  geometry  in  the  model.  However,  it  is 
probably  related  to  the  simplified  combustion  model.  The  breakup  of  the  jet  and  the  heat 
transfer  to  the  liquid  are  ignored  in  the  model.  Also,  the  pressure-dependent  burn  rate  used  is 
not  firmly  established,  especially  at  higher  pressures.  In  the  future,  variations  in  the  breakup 
and  combustion  models  will  be  examined.  Results  here  indicate  that  the  liquid  must  combust 
as  a  unit  in  a  reasonably  sized  volume  to  drive  the  first  radial  mode. 

8.  THE  155-mm  GUN 

Finally,  a  test  case  from  the  GE  first  generation  155-mm  gun  (Round  65)  is  modeled.  This 
is  a  7.2-liter  charge.  The  gun  fixture  and  this  particular  shot  are  discussed  in  detail  in  a 
previous  report  (Wren,  Coffee,  and  Morrison  1990).  The  lumped  parameter  model  was  run 
using  a  jet  breakup  model  developed  at  BRL  (Coffee  et  al.  1991).  The  piston  motion,  injection 
rate,  and  injected  droplet  size  from  the  lumped  parameter  model  were  used  as  input  into  the 
two-dimensional  code. 

Figure  36  compares  the  experimental  and  model  chamber  pressure.  The  model 
oscillations  start  a  little  later  and  die  out  much  sooner.  Analysis  indicates  that  the  oscillation 
die  out  because  the  injected  droplets  are  very  small.  The  liquid  burns  almost  instantaneously, 
and  there  is  no  liquid  accumulation  to  support  the  generation  of  oscillations. 

The  jet  breakup  model  used  was  developed  before  the  new  burn  rates  were  measured, 
and  the  slower  McBratney  rate  was  extrapolated  to  higher  pressures.  When  the  new  burn 
rate  is  used  at  high  pressures,  much  larger  droplets  can  be  injected  with  only  a  minor  effect 
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on  the  lumped  parameter  model.  So  the  lumped  parameter  model  was  modified  so  that  the 
smallest  droplets  injected  were  100  pm  in  diameter.  The  result  is  shown  in  Figure  37.  In  the 
model,  the  oscillations  still  start  a  little  late  and  die  out  too  rapidly.  However,  the  magnitude  of 
the  oscillations  is  a  very  good  match.  Figure  38  compares  the  experiment  and  model  at  the 
barrel  gauge,  B3,  which  is  uncovered  after  79.5  cm  of  projectile  travel.  The  model  gauge  is 
uncovered  a  little  too  soon,  and  the  oscillations  are  noticeably  smaller  than  in  the  experiment. 


In  Figure  39,  the  effect  of  refining  the  grid  is  considered.  The  top  curve  is  from  a  20  x  20 
grid  in  the  chamber  and  a  450  x  10  grid  in  the  tube.  The  middle  curve  is  from  a  40  x  40  grid 
in  the  chamber  and  a  900  x  20  grid  in  the  tube.  The  bottom  curve  is  from  a  60  x  60  grid  in 
the  chamber  and  a  1 ,000  x  30  grid  in  the  tube.  The  oscillations  are  very  similar  regardless  of 
the  grid  refinement.  It  is  possible  that,  if  the  grid  is  made  fine  enough,  the  oscillations  will 
diminish  in  magnitude.  However,  since  the  last  case  required  about  8  hr  of  run  time  on  the 
BRL  Cray  X-MP,  further  grid  refinement  is  impractical. 

The  oscillations  were  filtered  from  the  model  chamber  pressure,  and  the  filtered  pressure 
was  compared  with  the  lumped  parameter  model  (Figure  40).  The  agreement  is  very  good 
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Figure  38.  155-mm  Gun.  Shot  65. 
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Figure  40.  155-mm  Gun,  Shot  65.  Smoothed  Chamber  Pressure.  Lumped  Parameter  Model 


line).  Two-Dimensional  Model  (dot). 


until  near  the  peak  pressure,  when  the  two-dimensional  model  shows  a  little  lower  pressure. 
Since  the  lumped  parameter  model  has  a  simplified  algorithm  for  flow  into  the  gun  tube  and 
the  pressure  distribution  in  the  tube,  this  minor  disagreement  is  not  surprising.  However,  the 
projectile  velocity  almost  exactly  matches  the  experimental  data  (Figure  41 ). 

Figures  42-44  compare  the  Fourier  transforms  of  the  experiment  and  model  for  different 
time  windows.  The  agreement  is  actually  quite  reasonable.  Both  the  model  and  the 
experiment  show  a  broad  range  of  frequencies.  The  experiment  shows  higher  peaks  at  the 
higher  frequencies,  especially  at  later  times.  At  later  times,  the  model  tends  to  settle  in  at  a 
third  radial  mode  (at  approximately  1 0.5  kHz),  since  the  liquid  is  injected  about  a  third  of  the 
way  from  the  centerline  to  the  wall.  Figure  45  compares  the  transforms  for  the  barrel  gauge. 
The  agreement  is  good  for  the  lower  frequencies,  but  the  model  totally  fails  to  pick  up  the 
frequencies  from  1 0  to  20  kHz. 

One  of  the  problems  in  analyzing  the  experimental  data  has  been  the  unexpected  wide 
range  of  frequencies.  The  expected  acoustic  modes  of  the  chamber  are  fairly  small  (first 
radial  is  around  4  kHz).  In  liquid  propellant  rockets,  if  oscillations  occur,  they  generally  match 
the  lower  acoustic  modes  (Harrje  1972).  However,  the  operating  pressures  for  a  gun  are 
much  higher  than  for  a  rocket.  Large  amounts  of  liquid  have  to  be  injected  into  a  small 
volume.  The  energy  density  is  much  greater  than  for  a  rocket. 

One  proposed  mechanism  has  been  a  random  combustion  model  (Haberl  and  Pasko 
1 990).  That  is,  a  cloud  of  droplets  in  a  small  volume  ignites  and  bums  very  rapidly, 
generating  a  pressure  pulse.  Since  the  liquid  is  injected  rapidly,  a  dense  cloud  of  liquid  is 
expected  to  be  formed,  and  combustion  may  not  be  immediate.  Groups  of  droplets  could  be 
exposed  to  hot  gas  by  vortex  shedding,  or  combustion  in  a  volume  could  be  started  by  a 
previous  pressure  pulse.  For  a  previous  version  of  the  two-dimensional  code,  in  which  the 
two-phase  equations  were  not  implemented,  assuming  random  combustion  was  the  only 
mechanism  which  was  found  to  generate  high  frequencies  (Coffee  1 990). 

However,  with  the  present  version  of  the  model,  higher  frequencies  are  generated 
naturally,  without  the  necessity  of  putting  in  a  random  element.  The  liquid  is  injected  into  a 
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Figure  43.  155-mm  Gun,  Shot  65.  Chamber  Pressure.  Foi 
Experiment  (line^.  Model  (dot). 
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Figure  45.  155-mm  Gun.  Shot  65.  Gun  Tube  Pressure.  Fourier  Transforms.  14  to  15  ms. 

Experiment  (line).  Model  (dot). 

small  volume.  This  small  volume  combusts  rapidly  and  generates  a  high  pressure  region. 

The  pressure  is  kept  high  for  a  short  period  by  inertial  confinement,  and  a  pressure  pulse  is 
generated.  These  pressure  pulses  bounce  off  the  walls  and  interact  with  the  accumulated 
liquid. 

The  complicated  frequency  structure  is  probably  due  in  part  to  the  complicated  geometry 
of  the  chamber.  The  two-dimensional  model  was  also  run  without  the  gun  tube  model.  That 
is,  outflow  conditions  were  implemented  at  the  entrance  to  the  gun  tube,  with  the  pressure  just 
outside  the  chamber  obtained  from  the  lumped  parameter  code.  Figure  46  compares  the 
early  frequencies  from  the  model  with  and  without  the  gun  tube.  Without  the  gun  tube,  the 
frequencies  correspond  to  a  third  radial  mode  with  overtones.  With  the  gun  tube,  the  structure 
is  more  complicated.  The  tentative  conclusion  is  that  the  more  complicated  geometry  causes 
more  structure  in  the  frequencies.  In  the  actual  gun,  the  control  piston  extends  into  the 
chamber,  and  there  is  chambrage  from  the  chamber  to  the  tube. 

The  presence  of  oscillations  in  the  gun  tube  has  been  taken  as  evidence  that  liquid  must 
enter  the  tube  and  combust  there.  While  the  agreement  between  model  and  experiment  at 
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Figure  46.  155-mm  Gun.  Shot  65.  Chamber  Pressure.  Fourier  Transforms.  10  to  11  ms. 

Model  With  Tube  (line).  Model  Without  Tube  (dot). 

the  barrel  gauge  is  less  impressive,  the  model  does  show  that  pressure  oscillations  can  reach 
far  down  tube.  Figure  47  shows  the  liquid  volume  fraction  in  the  chamber  (volume  of  liquid 
over  total  control  volume).  The  graphs  shows  the  chamber  from  centerline  to  wall.  The  liquid 
combusts  well  before  the  entrance  to  the  gun  tube.  Figure  48  shows  the  corresponding 
pressure. 

Of  course,  the  model  greatly  simplifies  the  injection  process.  The  model  does  not  include 
the  actual  jet  breakup,  the  heat  transfer  to  the  liquid,  and  the  ignition  process.  Nevertheless, 
the  model  demonstrates  that  even  if  the  combustion  is  confined  to  the  chamber,  pressure 
waves  can  still  propagate  down  tube. 

9.  CONCLUSIONS 

A  two-dimensional  model  of  the  combustion  chamber/gun  tube  of  a  regenerative  liquid 
propellant  gun  has  been  developed.  The  code  has  been  extensively  checked  over  a  series  of 
test  problems.  The  code  has  been  compared  with  experimental  data  for  two  cases — a  Sandia 
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test  chamber  and  an  actual  155-mm  gun  firing.  While  agreement  between  model  and 
experiment  is  not  exact,  the  qualitative  features  of  the  experimental  data  are  reproduced. 

There  are  some  known  weaknesses  in  the  model.  The  geometry  is  restricted  to  simple 
cylinders.  The  liquid  and  gas  always  move  at  the  same  velocity  (infinite  drag).  The  actual 
injection  process  is  greatly  simplified  (that  is,  the  jet  breakup,  heat  transfer,  and  ignition 
process  is  not  included).  The  first  two  weaknesses  can  be  fixed.  However,  the  details  of  the 
injection  process  are  not  well  understood. 

It  has  long  been  a  puzzle  why  high  frequency  oscillations  are  generated  in  large  caliber 
liquid  propellant  guns.  The  present  model  indicates  this  is  related  to  the  rapid  injection  of 
large  amounts  of  liquid  propellant  into  a  small  volume.  Due  to  the  large  pressure  sensitivity  of 
the  burning  rate,  large  pressure  pulses  are  generated  just  by  inertial  confinement.  These 
pressure  pulses  bounce  off  the  walls  of  the  chamber  and  further  interact  with  the  unburned 
liquid.  In  liquid  propellant  rockets,  the  rate  of  injection  of  liquid  per  unit  volume  is  much 
smaller. 

There  are  several  possibilities  for  reducing  or  eliminating  the  oscillations.  If  the  burning 
rate  exponent  was  smaller,  oscillations  would  not  be  generated  because  the  combustion 
would  not  be  sensitive  to  pressure.  However,  no  known  reasonable  propellant  has  a  small 
enough  exponent.  Alternately,  the  burning  rate  could  be  increased,  so  that  the  liquid 
combusts  almost  as  soon  as  it  enters  the  chamber.  If  there  is  no  liquid  accumulation,  there  is 
no  way  to  generate  a  local  high  pressure  pulse.  This  can  be  done  mechanically,  by  breaking 
up  the  jet  into  small  droplets.  There  has  been  some  success  in  the  Sandia  test  fixture  using 
a  splitter  plate  just  inside  the  chamber  (Rychnovsky  1991).  Alternately,  the  propellant  could 
be  chemically  changed  to  bum  more  rapidly.  A  propellant  with  these  properties  is  about  to  be 
tested  (Klein,  Coffee,  and  Leveritt,  to  be  published).  Also,  there  are  devices  to  damp  out 
oscillations  used  in  liquid  propellant  rockets,  such  as  liners,  quarter-wave  chambers,  and 
Helmholtz  resonators.  The  model  presented  here  indicates  that  the  mechanisms  that  create 
the  oscillations  are  the  same  as  in  liquid  rockets,  and  the  same  damping  methods  should 
work.  Again,  work  with  the  Sandia  test  fixture  shows  some  success,  especially  with  liners. 
However,  the  moving  pistons  make  it  difficult  to  implement  these  designs  in  an  actual  gun. 
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LIST  OF  SYMBOLS 


b  -  covolume,  cm3/g 

-  liquid  sound  speed,  cm/s 
cG  -  gas  sound  speed,  cm/s 
c  -  mixture  sound  speed,  cm/s 
cp  -  specific  heat  at  constant  pressure,  J/g-K 

Cy  -  specific  heat  at  constant  volume,  J/g-K 

eL  -  chemical  energy  of  the  propellant,  J/g 
g0  -  conversion  constant  =  107g/MPa-cm-s2 
hL  -  enthalpy  of  the  liquid,  J/g 
hG  -  enthalpy  of  the  gas,  J/g 

k  -  turbulent  kinetic  energy,  cm2/s2 

K,  -  bulk  modulus  at  zero  pressure,  MPa 
Kj  -  derivative  of  the  bulk  modulus 

m  -  combustion  rate,  g/s 

Ml  -  liquid  mass,  g 

Mg  -  gas  mass,  g 

M  -  mixture  mass,  g 

p  -  pressure,  MPa 

Pk  -  rate  of  production  of  kinetic  energy,  cm2/s3 

R„  -  specific  gas  constant,  J/g-K 

S  -  surface  area,  cm2 

t  -  time,  s 

T  -  gas  temperature,  K 

u'  -  wall  friction  velocity,  cm/s 

v,  -  axial  velocity,  cm/s 

vy  -  radial  velocity,  cm/s 

VL  -  liquid  volume,  cm3 

VG  -  gas  volume,  cm3 

V  -  volume,  cm3 

V'  -  Lagrangian  volume,  cm3 

x  -  axial  distance,  cm 
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y  -  radial  distance,  cm 

e  -  turbulent  dissipation,  cm2/s3 

Eq  -  gas  volume  fraction 

Y  -  ratio  of  specific  heats 

pL  -  liquid  density,  g/cm3 

pQ  -  gas  density,  g/cm3 

p  -  mixture  density,  g/cms 

-  normal  stress,  g/cm-s2 

Xyy  -  normal  stress,  g/cm-s2 

-  shear  stress,  g/cm-s2 

p  -  dynamic  viscosity,  poise  =  g/cm-s 

p,  -  turbulent  viscosity,  poise  =  g/cm-s 

v  -  kinematic  viscosity,  cm2/s 
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